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SECTION  1 


INTRODUCTION 


Flow  fields  associated  with  shock  reflec¬ 
tions  often  times  contain  embedded  shear 
layers.  Typical  examples  are  slip  lines  from 
triple  points,  wall  boundary  layers,  and  wall 
jets  associated  with  double-Mach  reflec¬ 
tions.  Because  the  Reynolds  numbers 
associated  with  these  flows  are  very  large 
(typically  10%  per  centimeter  or  even 
larger),  such  shear  layers  are  essentially 
inviscid  tangential  velocity  discontinuities. 

It  has  been  known  for  more  than  100  years 
that  certain  types  of  subsonic  shear  layers 
are  inviscidly  unstable.  Lord  Rayleigh 
(1880)  proved  a  necessary  condition  and 
later  Tollmien  (1935)  proved  a  sufficient 
condition  for  the  unstable  growth  of 
perturbations  in  shear  layers.  This  is 
summarized  by  the  Rayleigh-Tollmien 
Inflection  Point  Theorem  which  states  that: 
solutions  of  the  inviscid  momentum  equa¬ 
tion  are  unstable  to  small  perturbations  if 
(and  only  if)  the  initial  shear  layer  profile 
contains  an  inflection  point  (Schlichting, 
1968). 

The  late-time  consequence  of  such  flow 
instabilities  is  that  perturbed  shear  layers 
eventually  roll  up  into  large-scale  vortex 
structures.  Today,  there  is  ample  exper¬ 
imental  evidence  that  shear  layers  evolve 
according  to  this  mechanism.  One  of  the 
most  famous  examples  is  the  experimental 
work  of  Brown  and  Roshko  (1974)  which 
demonstrates  that  shear  layers  roll  up  into 
organized  rotational  structures,'  and  that 
the  largest-scale  structures  are  essentially 
independent  of  the  flow  Reynolds  number 


(although  the  micro-scale  mixing  within  the 
large  structures  remains  Reynolds  number 
dependent).  About  the  same  time.  Winant 
and  Browand  (1974)  demonstrated  that  the 
growth  of  the  mixing  layer  is  governed  by 
the  pairing  mechanism  of  these  vortical 
structures.  Many  other  examples  of  co¬ 
herent  vortex  behavior  may  be  found  in 
M.  Van  Dyke’s  Album  of  Fluid  Motion 
(1982). 

Recently.  Monkewitz  and  Huerre  (1982) 
performed  a  linear  spatial  stability  analysis 
of  the  inviscid  momentum  equation,  and 
established  the  most  amplified  frequency 
(W|)  and  its  nth  subharmonics  (wn  =  W|/n) 
for  a  shear  layer  with  a  Tanh(y)  velocity 
profile.  Ho  and  Huang  (1982)  confirmed 
experimentally  the  effectiveness  of  these 
subharmonic  forcing  frequencies  in  modi¬ 
fying  the  dynamics  and  spreading  rate  of 
the  shear  layer.  Our  calculations  will  utilize 
both  results. 

Many  shear  layer  computations  involving 
discretizations  of  the  Navier-Stokes  equa¬ 
tions  are  available.  These  include  finite-dif¬ 
ference  methods  (e.g..  Corcos  and 
Sherman  |1984|,  Mclnville  et  al.  1 1985 1. 
and  Davis  and  Moore  [  1 985 1 ).  spectral 
methods  (Riley  and  Metcalfe.  1980).  and 
vortex  methods  (e.g..  Ashurst  ( 1 979 1 . 
Leonard  f  1 980 1 .  and  Ghoniem.  Chorin  and 
Oppenheim  1 1 982 1),  Oppenheim  |I986|). 
These  calculations  succeed  in  resolving 
small-scale  flowfield  structures  and  are  of 
high  quality.  However,  modeling  of  the 
viscous  terms  imposes  computational  inef¬ 
ficiencies  because  viscous  length  scales 


I.  Apparently  Brown  and  Roshko  were  not  the  first  to  observe  such  effects.  Michel  (1932) 
photographed  organized  rotational  structures  in  an  acoustically  perturbed  free  shear  layer. 
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must  be  resolved  on  the  mesh.  This  either 
severely  restricts  the  time  step  or  leads  to 
implicit  time  differencing  and  resultant 
numerical  diffusion.  The  random  vortex 
method  (Chorin  |1973|  and  |1986|)  avoids 
these  problems  at  the  expense  of  introduc¬ 
ing  a  statistical  error,  and  the  limitation  that 
the  flow  must  be  incompressible. 

Described  here  is  an  alternative,  large- 
Reynolds-nuniber  approach2  that  is  based 
on  the  nonsteady  solution  of  the  inviscid 
conservation  laws  of  gas  dynamics.  In 
contrast  with  the  vortex  methods,  there  is 
no  large-Mach-number  limit  with  this 
approach  and  the  baroclinic  generation  of 
vorticity  is  automatically  included.  In 
contrast  with  the  implicit  finite  difference 
solutions  of  the  Navier-Stokes  equations, 
the  present  approach  focuses  on  an 
accurate,  nondiffusive  evaluation  of  the 
convective  derivatives  which  seem  to 
control  the  dynamics  of  the  large-scale 
structures.  Hence,  molecular  diffusion 
effects  (i.e.,  diffusion  of  vorticity  and 
passive  scalars  across  stream-lines)  are 
neglected  during  the  time-scale  of  the 
calculation.  The  disadvantages  of  this 
inviscid  approach  are  that  calculations  of 
small-Mach-number  flows  are  expensive, 
and  shear  at  wall  boundaries  must  be 
programmed  into  the  calculation  as  an 
initial  condition. 

Numerical  results  were  obtained  by  means 
of  an  explicit  second-order  Godunov 
algorithm  (Colella  and  Glaz  [  1985a|)  which 


gives  nondiffusive  solutions  for  gas  dynam¬ 
ics.  This  algorithm  has  produced  accurate 
laminar  solutions  to  a  variety  of  blast  wave 
reflection  problems  (Colella  et  al.  |1985b| 
and  Glaz  et  al.  1 1 985 1 ).  Often  the  shear 
layers  embedded  in  these  flows  roll  up  into 
large-scale  vortex  structures  (Kuhl  et  al. 
1987).  One  of  the  questions  explored  here 
is:  How  accurate  are  these  nonsmooth  (i.e., 
fluctuating)  portions  of  the  flow,  as 
predicted  by  the  numerical  solutions  of  the 
inviscid  conservation  laws?  Illustrative 
calculations  are  presented  for  three  classes 
of  problems  typical  of  shock  reflection  flow 
fields:  free  shear  layers,  wall  boundary 
layers  and  a  wall  jet.  By  design,  they  are 
limited  to  nonsteady  calculations  of  two- 
dimensional  shear  layers  that  are  steady  in 
the  mean  or  time-averaged  sense.  The 
objectives  of  this  work  were  to  demon¬ 
strate  that  the  dynamic  evolution  of  these 
“steady”  shear  layers  can  be  predicted  by 
a  nonsteady  solution  of  the  inviscid  con¬ 
servation  law  of  gas  dynamics,  and  also  to 
evaluate  the  accuracy  of  such  solutions  by 
comparing  the  results  with  experimental 
data. 

The  numerical  formulation  of  the  calcula¬ 
tions  is  described  in  Section  2.  A  detailed 
discussion  of  the  results  for  the  free  shear 
layer,  wall  shear  layer  and  wall  jet 
calculations  are  presented  in  Sections  3.  4 
and  5,  respectively.  Conclusions  and  gen¬ 
eralizations  drawn  from  these  results  are 
offered  in  Section  6. 


2-  This  inviscid  approach  is  also  being  used  with  ihe  “Flux-Corrected  Transport"  codes  (Grinstein  el 
al.,  1986)  and  in  the  inviscid  vortex  dynamics  calculations  of  Inoue  (1985,  1987). 


2 


SECTION  2 


FORMULATION 


2.1  GOVERNING  EQUATIONS. 

The  time-dependent  flow  field  is  governed 
by  the  conservation  laws  for  mass,  momen¬ 
tum  and  total  energy.  Considered  here  is 
an  inviscid  gas  dynamics  approximation 
valid  for  large  Reynolds  number  flows; 
hence,  molecular  viscosity  and  heat  con¬ 
duction,  as  well  as  gravitational  forces  are 
neglected.  For  a  computational  volume  V 
with  boundary  jV,  the  equations  of  gas 
dynamics  may  be  written  in  strong  conser¬ 
vation  form  as: 

—  [  WdV 

r>t  J  y 

Wu-dA+l  aj  d/L  (1) 
M  J  dv  ~ 
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The  objective  of  a  numerical  algorithm  is 
to  evaluate  the  above  integrals  as  accurate¬ 
ly  as  possible. 

2.2  NUMERICAL  METHODS. 

The  computations  are  two-dimensional  and 
performed  in  Cartesian  coordinates.  A 
rectangular  mesh,  aligned  with  the  mean 
flow  (i.e.,  the  upstream  splitter  plate  or 
bottom  wall),  is  used.  Variable  mesh 
spacing  is  utilized  in  order  to  concentrate 
the  computational  effort  in  regions  of 
interest  and  minimize  any  adverse  effects 
of  the  artificial  outflow  boundary'  condi¬ 
tions.  All  computations  are  operator  split. 
The  temporal  approximation  of.  say.  an  x 
sweep  for  Eq.  (2)  may  be  written  in 
conservation  form: 


P  =  (Y-  1)  Pe 

In  the  above,  u  denotes  the  gas  velocity, 
while  p.  p  and  e  represent  the  gas  pressure, 
density  and  internal  energy,  respectively. 
For  convenience,  the  specific  heat  ratio  is 
assumed  here  to  be  constant:  y  =  1.4. 

Eq.  (1)  may  be  formally  integrated  over  a 
time  At  =  tn+l  -  tn  to  evaluate  new 
cell-averaged  values  of  the  conserved 
variables  W: 


(wv)  y 1  =  (wv) 


+ 


[cn+]/2  Cn+V2 

Fi*i/2,rFi-\/2.j 

\pn+\/2  pn+  1/2 

/+  1/1/2,;/-  1/2, y 


(3) 


where  the  cell  interface  fluxes.  F.  and 
pressure-work  terms,  P,  are  given  by: 
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c  n  +  1/2 

Fi^/2.j 


n  +  1/2  n  +  1/2 
i  +  1/2,7  +  1/2, y 


4  +  1/2,;  Ar 


n  +  1/2  n  +  1/2 
/  +1/2,y 55  2/  +i/2,; 


4+  1/2,7  A/ 


4+ 1/2,;  =  4(4+ 1/2,; ) 


These  interface  fluxes  and  pressures  are 
evaluated  by  the  second-order  Godunov 
technique  (Coiella  and  Glaz,  1985).  First, 
piecewise-linear  subgrid  interpolant  func¬ 
tions  are  used  to  define  the  environment 
everywhere  within  the  computational  cell  at 
the  initial  time  level.  tn.  Second,  the  slopes 
of  the  interpolants  are  limited  to  maintain 
monotonicity  of  the  solution.  Starting  at 
each  cell  boundary  (e.g.,  i+1/2,j)  at  time 
tn+l/2.  the  three  gas  dynamic  characteris¬ 
tics  are  traced  back  to  time  level  tn  where 
the  Riemann  variables  are  evaluated. 
Beginning  with  these  values,  the  Riemann 
equations  are  then  integrated  along  the 
characteristic  curves  to  evaluate  the  flow 
field  at  the  cell  interface  at  the  half  time 
level  tn+,/2.  Finally,  the  interface  fluxes 
and  pressures  based  on  this  solution  are 
used  in  the  conservative  difference  scheme 
of  Eq.  (3)  to  update  the  cell-averaged 
values  of  mass,  momentum  and  energy. 

2.3  TIME-STEP  CONTROL. 

The  above  represents  an  explicit  scheme 
for  solving  a  system  of  hyperbolic  conserva¬ 
tion  laws,  hence  the  time  step  is  limited  by 
the  Courant-Friedrichs-Lewy  condition: 


where  Ax  denotes  the  mesh  spacing,  and 
a  =  jyp/p  represents  the  local  sound 
speed.  The  Courant  number  used  in  these 
calculations  was  C  =  0.95. 

This  numerical  method  is  most  appropriate 
for  nonsteady  compressible  flows  of  gas 
dynamics.  It  can  also  be  used  to  provide 
inviscid  solutions  of  nonsteady  small 
Mach-number  flows,  but  one  must  either 
compromise  the  computational  efficiency 
or  violate  similitude  in  Mach  number.  For 
example,  the  Mach  numbers  typical  of 
many  shear  layer  experiments  are  on  the 
order  of  M  ~  0.01,  and  it  would  take 
approximately  one  hundred  time  steps  to 
convect  a  fluid  particle  through  one 
computational  cell  with  an  explicit  scheme. 
Note  that  there  is  nothing  theoretically 
wrong  with  this  approach.  In  fact,  it  can  be 
viewed  as  a  very  accurate  (albeit  expen¬ 
sive)  time  integration  scheme.  To  save 
computing  costs,  however,  the  sound  speed 
in  such  small  Mach-number  problems  was 
modified  so  that  flow  Mach  number  was 
typically  M  ~  0.2.  In  such  cases  the  density 
and  velocity  ratios  across  the  shear  layer 
were  preserved,  but  similitude  in  Mach 
number  was  violated.  For  low  speed  flows, 
however,  compressibility  or  Mach  number 
effects  are  generally  nor  important,  and  this 
Mach  number  compromise  was  found  to  be 
worth  the  improvement  in  computational 
efficiency. 

2.4  INTERFACE  DYNAMICS. 

To  clearly  illustrate  the  dynamics  of  shear 
layers,  it  was  necessary  to  accurately  track 
material  interfaces  (or  lines  in  2-D).  Such 
lines  were  discretized  by  a  series  of 
connected  points  ,t|<(t)  which  moved  with 
the  local  material  velocity,  i.e.: 

*  =  1.2 .  K  (?) 


At  =  CAx/(a  +  |//|) 


The  initial  spacing  between  points  was 
assumed  to  be  one-tenth  the  computational 
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cell  size.  The  location  of  each  material 
point  was  updated  according  to  a  finite 
difference  approximation  to  Eq.  (5), 
namely: 


place  the  boundaries  of  the  mesh  far  away 
from  the  layer  and  thereby  minimize  any 
boundary  influences  on  the  free  shear  layer 
development. 


n  + 1/2  n  n , 
x .  '  =X,+U  (x.(t  ))At/2 


(6) 


/7+t  /?  +  1/2  n,  (.n  +  1/2.,  At/~ 

■JU  =x.  '  +  u  (x.(t  '  ))Af/2 


where  un  represents  a  velocity  that  is 
linearly  interpolated  from  the  four  nearest 
computational  cells  surrounding  point  x^. 
Note  that  for  both  steps  in  Eq.  (6),  the 
velocity  field  is  based  on  time  tn  (because 
the  numerical  algorithm  is  operator  split); 
only  the  interpolation  point  changes:  i.e., 

x11^  versus  x^  +  When  the  distance  be¬ 
tween  neighboring  points 

exceeded  twice  the  initial  spacing,  a  new 
marker  was  inserted  at  their  midpoint  to 
maintain  numerical  resolution.  The  number 
of  marker  points  was  allowed  to  increase 
to  a  maximum  of  75,000  (which  would 
double  the  computational  time),  after 
which  no  further  points  would  be  inserted 
over  the  last  60  percent  of  the  grid.  Marker 
points  that  passed  outside  the  fine-grid 
region  were  deleted,  and  then  the  entire 
system  was  renumbered. 


2.5  INITIAL  CONDITIONS. 


The  calculations  were  performed  on  a 
two-dimensional  Cartesian  grid.  A  typical 
example  for  the  free  shear  layer  problems 
is  shown  in  Figure  1.  A  fine-zoned  region, 
covering  the  domain: 

0  «  x  «  Xp 
-  Yp  <  y  <  Yp 

was  centered  on  the  centerline  of  the  layer. 
An  expanding  grid  was  used  (above,  below, 
and  to  the  right  of  the  fine-grid  region)  to 


The  shear  layers  were  modeled  or  approxi¬ 
mated  by  a  Tanh(y)  velocity  profile.  This 
approach  not  only  allows  one  to  resolve  the 
shear  layer  on  the  computational  grid,  it 
also  removes  what  would  otherwise  be  an 
inviscid  velocity  discontinuity— -thereby  reg¬ 
ularizing  the  inviscid  problem.  In  effect, 
viscosity  is  allowed  to  act  on  the  flow  for 
a  long  enough  time  to  spread  the  shear  over 
a  finite  thickness  of  fluid;  this  solution  is 
then  mapped  onto  the  grid,  and  thereafter 
the  effects  of  molecular  viscosity  are 
eliminated. 

The  following  Tanh(y)  velocity  profile  was 
used  to  initialize  the  flow-field: 


u(x,y,o)  =fiy) 

(7) 

v(.t,  y,  o)  -  o 
where 

fiv)  =  Uml  1  +  A  Tanh  (  2 y/d  -  y0  )  |  (8) 

Um  =  «7,  +  U2)  /  2  (9) 

A  a  ((./)  -  f./2)  /  {U\  +  l>2)  (10) 

In  the  above,  U|  and  U2  denote  the  free 
stream  velocity  at  y  =  +00  respectively:  Um 
represents  the  mean  flow  velocity;  \  is  the 
shear  parameter:  yQ  denotes  the  centerline 
of  the  layer:  while  8  is  the  maximum-slope 
thickness  of  the  layer,  typically  2  cells.  In 
most  of  the  results  to  be  presented,  the 
lengths  have  been  normalized  by  8/2. 

In  some  free  shear  layer  problems,  the 
density  of  the  two  streams  is  different.  In 
these  cases  a  Tanh(y)  profile  was  also  used 
to  approximate  the  density  distribution 
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across  the  molecular  diffusion  region, 
namely: 

p  (.t,  y.  o)  =  h(y)  (11) 

where 

h(y)  =  Pm  1 1  +  2.  P  Tanh  (2y/<5  -  y0 )  j  (12) 
Pm  -  (  Pi  +  P2)  /  2  (13) 

Ap=(  Pj-P2)  /  (  P1  +  P2)  O4) 

In  the  above  pi  and  p2  represent  the  free 
stream  densities  at  y  =  +  00,  respectively; 
pm  denotes  the  mean  density  of  the  layer; 
and  \p  is  the  density  parameter, 

2.6  BOUNDARY  CONDITIONS. 

The  left  boundary  of  the  mesh  was  driven 
by  the  same  Tanh(y)  profiles  but  with  a 
sinusoidal  time  perturbation  on  the  stream- 
wise  velocity  component: 

11(0.  v,t)  =  u(.x,y,o)  *  g(t) 

(15) 

v(o.  y.  t)  =  o 

P(o,y,  t)  =  h(y)  (16) 

where 


N 

g(o  =  1  +  X  e" sin  ^  ( 1 1) 

n  =  1 

The  perturbation  frequencies  (wn)  and 
amplitudes  (€n)  were  taken  from  the  linear 
stability  analysis  for  Tanh(y)  shear  layers 
of  Monkewitz  and  Huerre  (1982): 

a)  1  =  0.219  and  u)n  =  o)\/n 

ei  =  0.01,  €2  =  0.75e], 

=  0.55ei  64  =  0.446^,  (18) 

and  e/7  =  1.7  for  n  >  5  (19) 

From  the  above,  that  the  maximum 
perturbation  amplitude  used  in  this  study 
was  one  percent  (e  1  =0.01).  Note  that  the 
above  disturbances  can  be  viewed  as 
modeling  the  spectrum  of  perturbations 
which  arise  naturally  in  the  flow  (e.g., 
molecular  fluctuations  which  are  amplified 
to  a  macroscopic  scale  by  means  of  linear 
stability  mechanisms). 

The  implementation  of  boundary  condi¬ 
tions  was  straightforward  for  the  calcula¬ 
tions  presented  here.  Dirichlet  conditions 
are  imposed  at  the  upstream  side  |  using 
eqs.  (15)  -  (17)  and,  for  the  free  shear  layer 
problems,  at  the  top  and  bottom  (using  free 
stream  values)  |.  The  downstream  side  is 
treated  with  standard  outflow  conditions  for 
all  variables.  Wall  boundary  layers  are 
treated  by  reflection,  thereby  allowing 
inviscid  slip. 
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SECTION  3 


FREE  SHEAR  LAYER  CALCULATIONS 


3.1  FORCED  SHEAR  LAYERS. 

The  purpose  of  these  calculations  was  to 
investigate  the  dynamic  response  of  a  shear 
layer  to  particular  perturbation  frequen¬ 
cies.  and  to  compare  the  results  with 
experimental  data.  To  that  end,  the  shear 
layer  parameters  were  chosen  to  corre¬ 
spond  to  the  forced  shear  layer  experiments 
of  Ho  and  Huang  (1982): 

U\  -  1.31  ,  U2  =  0.69  ,  Um  =  1.0 

A  =  0.31  and  r  =  U2/U\  =0.53 

Pi  =  p2  =  1  ,  A  p  =  o 

In  the  above,  the  velocity  components  have 
been  normalized  by  the  average  velocity  of 
the  two  streams,  and  the  density  by  the 
ambient  value.  The  layers  were  initialized 
with  uniform  thermodynamic  properties: 

?1  =  e2  -  44.63,  p\  *  p2  «  17.86, 

uj  =  a2  *  5.0,  y  »  1.4 

Note  that  by  this  choice  of  sound  speed,  the 
mean  flow  Mach  number  is  Mm  =  0.2.  The 
fine-grid  region  consisted  of  300  x-cells  by 
60  y-cells;  the  overall  mesh  covered  a 
domain  of  8600  by  3000.  The  initial  shear 
layer  width  was  2  cells  based  on  the 
maximum  slope  thickness,  or  say  80  =  6 
cells  based  on  the  99  percentile  change  in 
velocity.  A  one  percent  perturbation 
amplitude  was  used  for  each  frequency. 
The  dynamics  of  the  flow  was  elucidated 
by  tracking  a  material  line  that  was 
embedded  in  the  initial  center  line  (i.e., 
inflection  point)  of  the  shear  layer.  Each 


calculation  was  run  for  2000  cycles.  The 
final  interface  shapes  are  depicted  in 
Figure  2. 

In  the  Mode  I  calculation  (Fig.  2a).  only  the 
fundamental  frequency  W|  was  used  to 
perturb  the  layer.  Undulations  in  the  shear 
layer  are  first  noticeable  at  x  *  24.  while 
the  first  cusp  in  the  material  line  appears 
at  x  *  75.  At  larger  distances.  Figure  2 
shows  that  the  layer  rolled  up  into  discrete 
vortices  that  had  a  constant  wavelength  L| 
=  380  corresponding  to  the  fundamental 
frequency  u>i  (where  80  =  6).  The  width  of 
the  mixing  region  was  about  18„.  For  this 
Mode  I  response,  the  vortices  continued  to 
maintain  their  own  identity,  and  no  merging 
occurred  for  x  <  300.  There  was.  howev  er, 
some  “rebraiding  effects,"  where  material 
from  the  outer  bottom  edge  of  one  vortex 
structure  is  drawn  aft  and  down,  thus 
forming  the  lower  edge  of  the  trailing 
vortex  structure.  In  a  similar  fashion, 
material  from  the  top  of  the  trailing  vortex 
is  drawn  forward,  up  and  over  the  leading 
vortex  —  thereby  again  covering  the  braid 
region  between  vortices. 

In  the  Mode  II  calculation,  the  first  and 
second  frequencies  (W|  and  002)  were  used 
to  perturb  the  layer.  Initially,  vortices 
formed  at  the  fundamental  wavelength  L|. 
Soon,  however,  pairs  of  vortices  began 
interacting.  The  trailing  vortex  of  each  pair 
was  drawn  up  and  over  and  then  entrained 
into  the  leading  vortex  —  thus  forming  a 
larger-scale  structure.  Subsequently,  these 
structures  maintained  their  identity,  and  no 
further  merging  occurred  for  x  =s  300.  The 


Figure  2.  Numerical  simulation  of  the  forced  shear  layer  experiments  of  Ho  and  Huang 
(X  *  0.3 1  .r  =  0.53,  €  =  0.01).  Material  interface  plots:  (a)  Mode  I  response  (o», 
only):  (b)  Mode  II  response  (wi  +  w2);  (c)  Mode  Ul  response  (w(  +  W3);  (d)  Mode 
IV  response  (W|  +  w4). 
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merged  wavelength  was  twice  the  funda¬ 
mental  value  (L2  *  2L|  =  68c),  correspond¬ 
ing  to  the  second  frequency.  (02  =  u)|/2.  The 
final  width  of  the  mixing  region  was 
constant  at  about  2.580.  The  rollup  of  the 
material  line  in  the  merged  structures  was 
quite  smooth  and  basically  elliptical;  five 
complete  rotations  of  the  central  cusp  can 
be  observed.  The  start  of  rebraiding  effects 
are  visible  near  the  right  edge  of  the  grid. 

In  the  Mode  in  calculation,  the  first  and 
third  frequencies  (coi  and  (1)3)  were  used  to 
perturb  the  layer.  Just  as  before,  vortices 
initially  formed  at  the  fundamental  wave¬ 
length.  L|.  Soon,  however,  vortices  began 
interacting  in  groups  of  three,  and  a  more 
complex  merging  pattern  developed. 

The  third  (or  last)  vortex  was  drawn  up  and 
over  and  entrained  into  the  middle  vortex, 
forming  a  stronger  vortex  pair.  This  pair 
attracted  the  leading  vortex  which  was  then 
entrained  from  below.  The  wavelength  of 
the  final  merged  structures  was  three  times 
the  fundamental  value  (L3  =  3L[  *  980), 
corresponding  to  the  third  frequency,  0)3  = 
u)|/3.  The  final  width  of  the  mixing  region 
was  constant  at  about  58u-  Rebraiding  or 
further  merging  did  not  occur  for  x  <  300. 
The  rollup  of  the  material  line  in  the  large 
structures  was  quite  complicated  in  this 
case,  due  no  doubt  to  the  nonuniform 
distribution  of  vorticity  created  by  the 
merging  of  three  vortices. 

In  the  Mode  IV  calculation,  the  first  and 
fourth  frequencies  (uj|  and  0)4)  were  used 
to  perturb  the  layer.  Again,  vortices  formed 
at  the  fundamental  wavelength,  L|.  Soon, 
however,  vortices  began  interacting  in 
groups  of  four,  and  a  very  complex  merging 
pattern  developed.  First,  the  middle  two 
vortices  merged:  then  they  entrained  the 
leading  vortex  from  below;  finally,  the  last 


vortex  was  drawn  up  and  over  and 
entrained  into  the  main  structure.  The 
merged  wavelength  was  approximately  four 

times  the  fundamental  value  (I _ »  =  4L|  = 

1280).  The  computational  grid  was  too 
short  to  establish  the  asymptotic  vortex 
pattern  and  mixing  region  width.  As  in 
Mode  III,  the  rollup  of  the  material  line  in 
the  merged  structures  was  very  complicated 
due  to  the  nonuniform  distribution  of 
vorticity  created  by  the  merging  of  four 
vortices. 

The  above-described  dynamic  response 
duplicates  the  experimental  results  found 
by  Ho  and  Huang  (1982)  in  their  forced 
shear  layer  experiments  for  all  modes  that 
they  studied  Modes  1  through  IV). 

These  results  indicate  that  the  dynamic 
response  of  a  shear  layer  is  deterministically 
related  to  the  perturbation  frequencies.  If 
only  the  fundamental  frequency  (o>  1 )  and 
one  subharmonic  (u>n)  perturbations  are 
used,  the  vortices  merge  to  form  periodic 
large-scale  structures  with  a  constant  wave¬ 
length  Ln  =  nL|.  No  further  merging 
occurs,  and  the  width  of  the  mixing  region 
remains  constant  at  a  value  of  about 
2.5(n-l)80.  Hence,  continued  growth  of  the 
mixing  region  requires  the  addition  of 
smaller  and  smaller  frequency  (i.e.  longer 
and  longer  wavelength)  perturbations. 
Preliminary  calculations  supporting  these 
conclusions  were  reported  in  Chien  et  al. 
(1987). 

3.2  SPREADING  SHEAR  LAYER. 

A  variety  of  perturbation  frequencies  are 
present  in  most  actual  shear  layers,  and  the 
mixing  width  of  the  shear  layer  continues 
to  grow  with  increasing  distance.  This 
section  describes  a  numerical  simulation  of 
just  such  a  situation,  namely,  the  funda¬ 
mental  shear  layer  experiments  of  Brown 
and  Roshko  (1974).  The  particular  case 
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considered  was  that  of  a  5  m/s  helium  flow 
(state  1)  over  a  1.89  m/s  nitrogen  flow 
(state  2);  expressed  in  terms  of  the  shear 
layer  parameters,  these  become: 

U\  =  1.451,  U2  =  0.549,  Um*  1.00 
>1  =  0.451,  r  =  U2/U\  =  0.378 

Pi  *  1  ,  P2  =  7  ,  Ap  =-0.75 

In  the  above,  the  velocity  components  have 
been  nondimensionalized  by  the  average 
velocity  of  the  two  streams,  and  the  density 
by  p|.  By  design,  the  dynamic  pressure  of 
the  two  streams  are  equal. 

The  calculation  was  performed  on  a  grid 
with  a  fine-zoned  region  of  500  x-cells  by 
100  y-cells  (Ax  =  Ay  =  1 );  the  complete  grid 
covers  a  domain  of  8800  by  3000.  Again, 
the  initial  shear  layer  width  was  80  =  6  cells 
(based  on  the  99  percentile  change  in 
velocity).  The  flow  field  was  initialized  by 
a  Tanh(y)  profile  for  both  the  velocity  (Eq. 
8)  and  density  field  (Eq.  12),  and  a  uniform 
pressure:  pi  =  P2  =  17.86.  The  pressure  was 
selected  so  that  the  mean  flow  Mach 
number  was  Mm  =  0.4  (while  Mt  =  M2  = 
0.29).  Only  the  velocity  field  was  perturbed 
on  the  left  boundary  (according  to  Eq.  15). 
For  the  first  thousand  cycles  only  the 
fundamental  perturbation  frequency  wj 
was  used  so  that  the  mixing  started 
gradually;  thereafter  (cycles  1,000  to 
10.000)  the  first  ten  modes  were  used  to 
model  the  natural  growth  of  a  spectrum  of 
small  perturbations  present  in  most  experi¬ 
ments.  A  material  line  (initially  on  the 
center  line  of  the  shear  layer)  was  tracked 
to  follow  the  distortion  of  the  interface. 

The  dynamic  evolution  of  the  shear  layer 
is  depicted  in  Figure  3  by  means  of  material 
interface  plots  at  various  times.  The  first 
cusp  or  kink  in  the  material  line  appears 


at  x  s  65.  Clockwise-rotating  vortex 
structures  form  initially  at  the  fundamental 
wavelength  Lj;  but  these  rapidly  merge  to 
form  larger-scale  structures.  In  this  way, 
the  mixing  width  continues  to  grow  with 
increasing  streamwise  distance. 

The  response  of  the  shear  layer  was 
approximately  periodic.  Figure  3  depicts 
the  evolution  of  the  interface  over  one  such 
period  from  frame  a  (t  =  4035)  to  frame 
s  (t  =  501 3) |;  frame  t  shows  that  this 
interface  shape  reappears  again  at  t  =  6018, 
thus  the  period  of  the  overall  flow  is  At  m 
1000.  The  following  is  a  detailed  descrip¬ 
tion  of  the  development  of  the  large-scale 
structures  depicted  in  frame  s  (t  =  5013). 
The  individual  vortex  structures  are  labeled 
sequentially  with  cardinal  numbers.  Con¬ 
sidered  first  is  the  formation  of  the 
largest-scale  structure,  of  frame  s,  structure 
6-12.  Vortices  6,  7  and  8  merge  (frames  d 
through  g)  in  a  Mode  III  response  to  form 
structure  6-8.  Vortices  9  and  10  pair  and 
merge  (frames  g  through  j)  in  a  Mode  II 
response;  vortices  11  and  12  also  pair  and 
merge  (frames  h  through  j).  Finally, 
structures  6-8  and  9,  10  and  11.12  interact 
and  merge  to  form  the  largest-scale 
structure  6-12  shown  in  frame  s. 

In  the  meantime,  structure  13-18  has  been 
developing.  Vortices  13,  14  and  15  merge 
(frames  i  through  I)  in  a  Mode  III  response; 
vortices  16,  17  and  18  also  merge  (frames 
j  through  I)  in  a  Mode  III  response.  Then 
structure  13-15  pairs  and  merges  with 
structure  16-18  (frames  n  through  p)  in  a 
Mode  VI  response,  resulting  in  the 
large-scale  structure  13-18. 

At  the  same  time,  structure  19-24  has  been 
developing.  Vortices  19,  20  and  21  merge 
(frames  k  through  m)  in  a  Mode  III  re¬ 
sponse;  vortices  22,  23  and  24  also  merge 
(frames  m  through  o)  in  a  Mode  ill 
response.  Finally,  structure  19-21  pairs 
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Figure  3.  Numerical  simulation  of  the  spreading  shear  layer  experiment  of  Brown  and 
Roshko  (\  =  0.451,  r  =  0.378,  p2/pi  =  7)  (Continued). 
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with  structure  22-24  (frames  r  and  s)  in  a 
Mode  VI  response,  resulting  in  structure 
19-24  as  shown  in  frame  s. 

In  summary,  the  large-scale  structures  are 
created  by  Mode  II,  III  and  VI  merging 
processes.  Interface  shapes  periodically 
reappear  every  At  »  1000. 

The  material  interface  in  Figure  3  is 
actually  represented  by  unconnected  dots 
or  tracer  points.  Initially  the  interface 
contained  5,000  such  points  (ten  points  per 
Ax).  The  number  of  tracers  rapidly 
increased  during  the  calculation  to  a 
maximum  of  75,000  points  due  to  the  large 
strain  in  the  braid  regions  between 
structures.  The  length  of  the  interface 
increased  by  at  least  a  factor  of  15  due  to 
the  continued  vortex-merging  process. 

The  initial  material  interface  fluid  seems  to 
cover  a  significant  fraction  of  the  cross-sec¬ 
tional  area  of  the  large  structures.  This  is 
a  consequence  of  mixing  by  an  inviscid, 
nondiffusive  folding  of  the  material  inter¬ 
face.  Hence,  regions  of  essentially  pure 
helium  are  interleaved  with  regions  of  pure 
nitrogen.  This  effect  is  visible  in  exper¬ 
imental  photographs  to  be  shown  next. 

Figure  4(a)  depicts  the  calculated  material 
interface  at  t  *  5013.  For  comparison, 
Figure  4(b)  presents  a  shadowgraph  picture 
of  the  helium-nitrogen  interface  recorded 
during  the  shear  layer  experiments  (Brown 
and  Roshko.  1974;  Fig.  3d).  Similarities 
between  the  calculated  and  experimental 
interface  are  remarkable.  The  calculated 
spreading  rate  varied  between  18  and  22 
percent,  which  is  in  excellent  agreement 
with  visual  spreading  rate  of  8'Vis  =  21 
percent  observed  in  the  experiment.  The 
shape  and  wavelength  of  the  large-scale 
structures  are  quite  similar.  The  fine  white 


and  black  lines  on  the  large  structures  of 
the  experiment  may  now  be  understood,  in 
light  of  Fig.  4a,  as  folded  or  interleaved 
regions  of  pure  helium  or  pure  nitrogen. 


Figure  5  depicts  the  flow  field  near  the 
center  line  (i.e.,  y  *  0.5)  at  t  =  5013.  Large 
fluctuations  are  evident  in  the  velocity, 
density  and  dynamic  pressure.  In  fact, 
densities  corresponding  to  pure  helium  and 
to  pure  nitrogen  are  seen  near  the  center 
line;  similarly,  the  velocity  fluctuates 
between  values  corresponding  to  approxi¬ 
mately  U|  and  to  U2.  Overpressure  fluc¬ 
tuations  are  small  (between  +1.5  to  -4.0 
percent  of  ambient),  which  is  quite 
comforting  in  light  of  the  Mach  number 
compromise  (0.29  <  M  <  0.4)  used  for  this 
calculation. 


Also  included  in  Figure  4  are  the 
corresponding  density,  vorticity  and  over¬ 
pressure  contours.  The  shape  of  the  large 
structures  as  deduced  from  the  density 
contours  agrees  with  the  material  interface 
shape,  indicating  little  numerical  diffusion. 
It  is  interesting  to  note  that  vorticity  is 
created  near  the  braid  regions  and 
entrained  into  the  structures.  Approximate¬ 
ly  circular  low  pressure  regions  are  found 
in  each  large  structure,  and  this  generates 
a  radial  pressure  gradient.  As  the  he¬ 
lium-nitrogen  interface  is  entrained 
obliquely  through  this  pressure  gradient, 
vorticity  is  generated  by  the  baroclinic 
| V  (1/p)  x  V  p|  mechanism.  The  magnitude 
of  the  baroclinicly-generated  vorticity 
reaches  about  minus  one  half  of  the  inflow 
vorticity,  hence  it  must  have  an  effect  on 
the  flow.  Indeed,  close  inspection  of  the 
material  interface  reveals  some  small- 
er-scale  vortices  rotating  in  the  counter¬ 
clockwise  direction.  This  causes  the 
merging  patterns  to  be  considerably  more 
complex  and  less  regular  than  in  Figure  2. 
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Figure  4.  Flowfield  of  the  spreading  shear  layer  calculation  at  t  =  5013:  (a)  material 
interface:  (b)  shadowgraph  from  the  Brown  and  Roshko  experiment  (Figure  3): 
(c)  density  contours:  (d)  vorticity  contours  (solid  lines  denote  negative  values 
and  dashed  lines  correspond  to  positive  values  that  are  baroclinicly  generated): 
(e)  overpressure  contours  (solid  lines  denote  positive  values). 
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Figure  5.  Flow  field  along  the  center  line  (y  =  0.5)  of  the  spreading  shear  layer  calculation 
at  t  *  5013:  (a)  density;  (b)  streamwise  velocity;  (c)  overpressure:  (d)  dynamic 
pressure. 
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Monitoring  stations  located  at  x  =  100,  200, 
and  400  were  used  to  store  the  flow  field 
time  histories.  These  were  then  integrated 
in  time  over  the  last  2000  cycles  of  the 
calculation  (about  6  periods  of  the 
large-scale  structures  at  x  =  400)  to 
establish  the  mean  shear  layer  profiles.  The 
mean  velocity  and  density  profiles  are 
depicted  in  Figure  6.  The  three  stations  are 
compared  by  means  of  the  scaling  variable: 

nsL=y  /  (■*-*>)  (2°) 

where  x()  denotes  the  effective  origin  of  the 
breakdown  of  the  shear  layer  (x0  =  25, 
according  to  Fig.  5).  The  profiles  seem  to 
collapse  reasonably  well  with  this  scaling. 
The  shear  layer  is  resolved  across  approxi¬ 
mately  10,  25  and  50  computational  cells 
at  the  first,  second  and  third  stations, 
respectively. 

The  shaded  curves  in  that  figure  denote  the 
mean  profiles  as  measured  by  Brown  and 
Roshko  (1974).  The  calculated  density 
profiles  are  in  excellent  agreement  with  the 
experimental  profiles  —  even  to  such 
detailed  features  as  multiple  inflection 
points  and  plateau  regions  (e.g.,  near  t)  = 
0.05).  Agreement  between  the  calculated 
and  measured  mean  velocity  profiles  is  also 
good,  considering  the  accuracy  of  evaluat¬ 
ing  the  velocities  from  pitot  tube  measure¬ 
ments  for  flows  with  large  density 
variations. 

The  mean  transverse  velocities  lie  in  the 
range: 

r  0.016  (for  >  0.05) 

v/U |  s  '  0  (for//5L  =  o) 


pressure  profile  was  constant  and  equal  to 
the  ambient  value. 

Figure  7  presents  the  fluctuating  flow 
profiles  corresponding  to  the  mean  flow 
depicted  in  Figure  6.  These  were  calculated 
as  root-mean-square  time  averages  about 
the  mean  values  using  the  last  2000  cycles. 
There,  u  ’  and  v  ’  represent  the  streamwise 
and  transverse  velocity  fluctuations,  u’v' 
denotes  the  fluctuating  shear  stress,  while 
p’  and  p’  represent  the  density  and 
pressure  fluctuations,  respectively.  The 
fluctuating  flow  profiles  are  reasonably 
smooth,  except  for  u’v'  which  probably 
requires  a  longer  time  average  to  achieve 
convergence.  The  profiles  are  remarkably 
similar,  considering  that  the  rotational 
structures  at  x  =  400  are  much  more 
complex  than  at  x  *  100.  The  shaded  region 
in  Figure  7d  denotes  the  fluctuating  density 
profile  measured  by  Konrad  (1977)  for  the 
same  shear  layer  that  was  studied  by  Brown 
and  Roshko.  The  calculated  density  profiles 
are  in  excellent  agreement  with  this  data 
—  even  the  detailed  features  such  as  the 
peak  values  and  the  plateau  region  near 
t)sl  =  0.07  agree  quite  well. 

The  dashed  curves  in  Figure  7  represent  the 
r.m.s.  velocity  fluctuations  as  measured  by 
Oster  and  Wygnanski  (1982:  their  Figs.  6a. 
7b  and  7c)  for  a  constant  density  shear 
layer  with  a  different  velocity  ratio.  Our 
calculated  profiles  are  only  qualitatively 
similar  to  their  profiles.  For  example,  they 
find  only  single-peaked  profiles,  while  our 
calculated  profiles  of  u'  and  u'v'  have  a 
central  peak  and  side  lobes,  the  latter  being 
caused  by  the  large  density  in  our  problem. 


v- 0.008  (for  <- 0.05) 

indicating  that  the  mixing  creates  an 
outward  displacement  of  the  streamlines  on 
either  side  of  the  shear  layer.  The  mean 


The  peak  values  of  the  fluctuating  quanti¬ 
ties  are  compared  with  other  shear  layer 
studies  in  Table  1.  The  present  results 
agree  quite  well  with  the  peak  values 
obtained  from  the  2D,  inviscid.  vortex 
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Figure  6.  Mean-flow  profiles  of  the  spreading  shear  layer  calculation:  (a)  streamwise 
velocity  profiles;  (b)  density  profiles.  Symbols  □.  O  and  A  denote  stations  at 
x  =  100,  200  and  400.  respectively.  Shaded  area  denotes  the  experimental  data 
band  from  Brown  and  Roshko  (1974). 
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Figure  7.  Time-averaged  fluctuating  flow  profiles  from  the  spreading  shear  layer  calcu¬ 
lations:  (a)  streamwise  velocity;  (b)  transverse  velocity;  (c )  shear  stress:  (cl) 
density  (shaded  region  corresponds  to  the  experimental  data  band  of  Konrad 
1 1 977 1 );  (e)  pressure.  Symbols  □.  O  and  A  denote  stations  at  x  =  100.  200  and 
400.  respectively.  Dashed  curves  denote  the  data  of  Oster  and  Wygnanski  ( 1982) 
for  an  unforced  shear  layer  (their  Region  1). 
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Table  1.  Comparison  of  Peak  Fluctuating  Quantities. 


Case 

Conditions 

Peak  Values  (10  2) 

FREE  SHEAR  LAYERS 

X 

r 

xp 

d  /U  i 

v'  /Uj 

w'  /U  i 

-if  v' /U,2 

P'  /P2 

P'  /pi 

Vortex  Dynamics  Calc 

0.25 

0.6 

0 

13 

9.3 

_ 

0.32-0.3  8 

_ 

(lnoue,1985) 

Oster  and 

0.25 

0.6 

0 

7.2 

6.1 

5.8 

0.21 

— 

— 

Wygnanski  (1982) 

0.45 

0.38 

0 

11.2* 

9.5a 

9.0a 

0.50a 

— 

— 

Present  Calculations 

0.45 

0.38  - 

-0.75 

x  =  100 

9.4 

11.2 

— 

0.54 

29 

1.3 

x  =  200 

15.9 

15.4 

— 

0.25 

35 

1.5 

x  =  400 

14.5 

18.8 

— 

0.69 

35 

1.7 

WALL  SHEAR  LAYERS 

X 

r 

^e 

U  /U<* 

V'/Uoo 

w'/Uoo 

-u' v'/U» 

P’/p» 

Shock  Tube  Case 

x  =  400 

-0.39 

2.27 

0.2C 

19 

10 

— 

0.13 

11 

8 

x  =  600 

26 

12 

— 

0.16 

13 

9 

x  =  800 

29 

14 

— 

0.12 

12 

10 

x  =  950 

31 

15 

— 

0.059 

10 

11 

Owen  et  al.  (1975) 

— 

— 

— 

— 

10 

— 

Tripped  Flat  Plate 

Case 

1 

0 

0 

x  =  400 

18.5 

7.26 

— 

0.0531 

0.169 

0.212 

x  =  600 

18.1 

7.36 

— 

0.0259 

0.156 

0.191 

x  =  800 

17.6 

7.25 

— 

0.0226 

0.144 

0. 1 75 

x  =  950 

17.1 

7.15 

— 

0.0186 

0.135 

0.162 

Klebanoff  (1955) 

8.75 

4.0 

6.4 

0.1 38b 

— 

— 

WALL  JET 

X 

r 

xp 

u'/Uj 

v'/Uj 

Present  Calculations 

-1 

OO 

0 

x  =  85  (5H) 

23 

28 

— 

1.85 

0.75 

1.0 

x  -  170  ( 10H) 

25-28 

33 

— 

4.5 

0.83 

1.1 

Bajura  and 

-1 

90 

0 

Catalano  (1975) 

x  =  10H 

0.7 

— 

— 

— 

— 

— 

x  -  20H 

2.0 

— 

— 

— 

“Scaled  values,  based  on  AL'/U\  =  0.622. 
hl'valualcd  from  llq.  38 
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dynamics  calculations  of  Inoue  (1985).  The 
present  peak  values  are  1 .5  to  2  times  larger 
than  the  experimental  data  of  Oster  and 
Wygnanski  (1982).  Of  course,  there  are 
differences  in  the  problems  considered 
(e.g..  different  X  and  Xp),  and  the 
experimental  data  certainly  includes  the 
effects  of  viscous  diffusion  and  dissipation, 
but  the  latter  effects  should  be  small 
because  of  the  large  Reynolds  number  of 
the  flow.  Probably  the  main  reason  why  the 
calculated  peaks  (especially  u '  and  v ' )  are 
larger  than  the  data  is  that  the  calculations 
are  two  dimensional  while  the  experimental 
data  contained  three  dimensional  fluctua¬ 
tions  —  that  is,  some  of  the  fluctuating 
kinetic  energy  is  stored  in  the  third 
dimension. 

3.3  SUMMARY. 

These  results  indicate  that  the  dynamic 
response  of  a  Ji  .urbed  shear  layer  is 
basically  independent  of  the  initial  shear 
layer  profile,  and  independent  of  the 
perturbing  frequency  as  long  as  a  spectrum 
of  frequencies  are  available  to  the  flow. 
Small  perturbations  cause  the  shear  layer 
tc  roll  up  into  periodic  vortices.  Subhar¬ 


monic  perturbations  cause  these  vortices  to 
merge  —  forming  successively  larger-scale 
rotational  structures.  The  dynamics  of  the 
large-scale  structures  determines  the  mix¬ 
ing  of  the  two  streams  and  controls  the 
fluctuations.  The  present  inviscid  numeri¬ 
cal  simulations  agree  rather  well  with 
experimental  data  —  not  only  for  such 
macroscopic  features  like  the  spreading 
rate  (8'vjS  *  0.21)  and  the  shape  and 
wavelength  of  the  large-scale  structures, 
but  also  in  mean  flow  profiles  such  as 
density  and  velocity.  Even  the  time- 
averaged  fluctuating  flow  profiles  such  as 
density  agree  with  the  data.  Hence,  we 
conclude  that  the  dynamic  evolution  of  the 
large-scale  structures  is  dominated  by 
inviscid  effects  that  are  well  approximated 
by  nonsteady  solutions  of  the  conservation 
laws  of  gas  dynamics. 

The  present  large-Reynolds-number  ap¬ 
proach  (i.e.,  the  direct  solution  of  the 
nonsteady  conservation  laws  of  gas  dynam¬ 
ics)  offers  then  one  technique  forextending 
the  numerical  simulations  of  fluctuating 
shear  layers  to  the  compressible  flow 
regime.  The  next  section  presents  an 
example  of  just  such  a  case. 
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SECTION  4 


WALL  SHEAR  LAYER  CALCULATIONS 


4.1  SHOCK  TUBE  EXPERIMENT. 

Considered  next  is  the  evolution  of  a  wall 
shear  layer  (i.e..  a  boundary  layer)  behind 
a  shock  wave  in  a  shock  tube.  The  key 
physical  aspects  of  the  problem  are 
depicted  in  Figure  8.  This  is  a  tracing  of 
a  shadow  photograph  provided  by 
Dr.  H.  Reichenbach  (Ernst  Mach  Institut, 
FRG).  A  planar,  constant  velocity  shock 
(Mach  number  Ms  =  1.7)  is  shown  prop¬ 
agating  to  the  left.  The  test  section  floor 
consisted  of  a  smooth  steel  plate  (rough¬ 
ness  height  less  than  I  pm),  while  the  roof 
consisted  of  a  smooth,  porous  plate 
(ceramic  Filtrokelit;  roughness  height  less 
than  30  pm).  Periodic  acoustic  waves 
emanate  from  the  point  where  the  shock 
front  intersects  the  floor  or  roof.  These 
waves  have  a  constant  spacing  of  about 
1  mm.  Scanning-electron-microscope  pho¬ 
tographs  of  the  surface  show  that  even  such 
nominally  smooth  surfaces  are  still  rough 
on  length  scales  which  are  large  compared 
to  the  shock  front  thickness.  Hence,  a  shock 
wave  cannot  propagate  with  uniform 
velocity  very  near  the  surface.  It  reflects 
from  the  front  face  of  each  roughness 
element,  thereby  enhancing  its  strength, 
and  as  it  propagates  down  the  backside  of 
each  roughness  element,  it  loses  strength. 
Therefore,  near  the  surface  the  shock  front 
velocity  must  oscillate,  even  though  it 
propagates  with  a  constant  velocity  at  large 
distances  from  the  surface. 

Figure  8  also  depicts  periodic  density 
structures  in  the  roof  and  floor  boundary 
layer.  (The  wall  boundary  layer  gas  is 
cooler  and  denser  than  the  free-stream 
flow,  and  therefore,  detectable  by  schlieren 


or  shadowgraph  techniques).  These  orga¬ 
nized  structures  are  inclined  to  the  flow  at 
an  angle  of  approximately  45°.  Their  wave¬ 
length  is  regular  and  constant  at  about  1 
mm,  which  is  similar  to  the  acoustic  wave 
spacing  and,  of  course,  much  longer  than 
any  wall  roughness  wavelengths.  Similar 
organized  structures  have  been  reported  by 
Spina  and  Smits  (1987)  for  turbulent 
boundary  layers  in  steady,  compressible 
flow  (while  Prandtl  [  1 933 1 .  Bergh  1 1 958 1 
and  Gad-el-Hak  et  al.  1 1 984 1  have  reported 
the  existence  of  organized  structures  in 
steady,  incompressible  boundary  layers). 
Hence,  we  postulate  that  the  oscillations  in 
the  shock  front  velocity  trip  the  wall 
boundary  layer,  which  rolls  up  into  periodic 
vortex  structures  and  thereby  generates 
acoustic  waves  that  radiate  away  from  the 
foot  of  the  shock  front. 

The  auxiliary  scale  in  Figure  8  depicts  the 
Reynolds  number  based  on  distance  behind 
the  shock  front.  The  flow  reaches  a  critical 
Reynolds  number  of  3  x  105  at  a  distance 
of  only  15  mm  behind  the  shock,  hence, 
essentially  the  entire  boundary  layer  may 
be  characterized  as  large  Reynolds  number, 
turbulent  flow. 

The  flow  conditions  across  the  Ms  =  1.70 
shock  front  are  listed  in  Table  2.  In 
stationary  coordinates,  states  1  and  2 
denote  conditions  ahead  and  behind  the 
shock  front,  respectively. 

4.2  FORMULATION. 

A  numerical  simulation  of  this  shock  tube 
experiment  was  performed  in  the  same 
spirit  as  the  simulation  of  the  Brown  and 
Roshko  shear  layer  case.  Namely,  we  test 
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Figure  8.  Tracing  of  a  shadow-schlieren  photograph  of  the  waves  induced  by  the  boundary 
layer  behind  a  Mach  M,  =  1.7  shock  wave.  Acoustic  waves  radiate  from  the  foot 
of  the  shock.  Density  structures  are  formed  in  the  wall  boundary  layer.  (Courtesy 
of  Dr.  H.  Reichenbach,  Ernst  Mach  Institute  FRG.) 
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Table  2.  Shock  Conditions  (Ms  =  1.70,  7  =  1,4). 


~  1  ••  1 

Stationary 

Shock-Fixed 

Coordinates 

Coordinates 

Wall 

Edge 

Flow  Field 

State 

1  State  2 

(  y=0  )  (  y®*) 

p/pi 

1 

3.2 

~3.2 

•3.2 

p/pi 

1 

2.2 

3.2 

2.2 

e/ei 

1 

1.46 

1.0 

1.46 

a/ai 

1 

1.21 

1.0 

1.21 

u/ai 

0 

0.94 

1.70 

0.77 

u/a  (local) 

0 

0.77 

1.70 

0.64 

the  postulate  that  the  dynamics  of  this  wall 
boundary  layer  flow  is  dominated  by  two- 
dimensional  inviscid  flow  effects  by 
comparing  the  calculated  results  with 
boundary  layer  data. 

The  calculation  was  performed  in 
shock-fixed  coordinates  with  the  shock 
located  at  x  =  o.  similar  to  those  used  by 
Mirels  (1956).  In  these  coordinates,  the  gas 
velocity  varies  between  the  value  ue  at  the 
edge  of  the  boundary  layer  (y  =  ») 

lie/o  1  a  (u'5-»2)/wi  =0.77 

to  a  maximum  value  uw  on  the  wall  (y  = 
o) 

=  1.70 


nates.  The  temperature  varies  between  the 
edge  value  of  430  K,  to  the  wall  value  of 
294  K,  with  the  latter  corresponding  to  an 
ambient  temperature  wall  boundary  condi¬ 
tion.  Other  flow  variables  are  listed  in 
Table  2. 

The  calculation  was  performed  on  a  grid 
with  a  fine-zoned  region  consisting  of  500 
x-cells  (Ax  =  2)  by  60  y-cells  (Ay  =  1).  as 
shown  schematically  in  Figure  9.  Expand¬ 
ing  cells  were  used  above  and  to  the  right 
of  the  fine-zoned  region  to  cover  a  total 
domain  of  17600  by  3000. 

The  laminar  boundary  layer  behind  the 
shock  was  approximated  by  a  Tanh(y) 
profile  for  the  streamwise  velocity  and 
internal  energy: 


where  the  gas  moves  with  the  shock 
velocity,  the  latter  being  precisely  the 
no-slip  condition  in  shock-fixed  coordi¬ 


u(.x,y,o)/ul  =  Ms  fly) 
v(.t,y,  0)  =  0 


(21) 


(a)  Stationary  coordinates 


(b)  Shock  •  fixed  coordinates 


where 


e(.r,y,  o)/e\ 

=  1 .22  +  0.244  Tanh(0.5y  -  1 .5)  (22) 

where 
f  (  y  ) 

=  0.741  -  0.286  Tanh{0.5y  - 1.5)  (23) 

For  convenience,  the  flow  variables  were 
nondimensionalized  by  the  ambient  condi¬ 
tions  denoted  by  subscript  1.  By  this 
formulation,  u/ai  =  Ms  =  1.70  on  the  wall, 
and  u/a|  =  0.455  Ms  =  0.77  in  the  free- 
stream,  thereby  introducing  the  correct 
amount  of  circulation  into  the  flow. 
Similarly,  e/ei  =  1  on  the  wall  and  e/e\  = 
1.46  in  the  free-stream,  thereby  satisfying 
the  correct  temperature  boundary  condi¬ 
tions.  The  pressure  field  corresponded  to 
the  uniform  state  behind  the  shock: 

p(x.y.o)/P]  =  1  +  l)2y/(y+  1)  =  3.2 

(24) 

while  the  density  profile  was  related  to  the 
pressure  and  internal  energy  profile  by 
means  of  the  equation  of  state: 

P  -  P/ e(y  ~  1 )  (25) 

This  gives  a  density  of  p/p  j  =  3.2  on  the  wall 
and  p/pl  =  2.2  in  the  free-stream. 

The  above  conditions  were  used  to  initialize 
the  flow  field.  The  left  boundary  of  the 
computational  mesh  was  then  driven  with 
a  sinusoidal  perturbation  on  the  streamwise 
velocity  profile: 

ti{o.y.t)/a\  -  A/(y,  t) 

(26) 

v(o,  y,  t)  =  0 
e(o.  y,  t )  =  e(.x,  v,  o) 

(27) 

p(o.y,0  =p/e(y-  1) 


10 

A/(y,  0  =  Ay)  1  =  <p(y)  2  €"  Sina>nt 

n  =  1 

(28) 

fl  -y/60  for  y  «  60 

<t>(y)  =  (29) 

0  for  y  >  d0 

By  this  formulation,  the  maximum  pertur¬ 
bation  occurs  on  the  wall,  and  is  zero 
outside  the  shear  layer  (y  >  8„  =  7.65)  — 
thus  modeling  the  physical  problem  as 
described  in  the  previous  section.  The 
perturbation  frequencies  a)n  were  assumed 
to  be  the  same  as  for  the  free  shear  layei 
case  (Eqs.  18),  although  it  was  recognized 
that  there  can  be  some  (as  yet  undeter¬ 
mined)  change  in  frequencies  due  to  the 
influence  of  the  wall.  A  maximum  pertur¬ 
bation  amplitude  of  6|  *  0.01  (with  €n  from 
Eqs.  19),  was  used,  corresponding  to  a  + 
1  percent  variation  in  shock  velocity  near 
the  wall.  The  pressure  was  coupled  to  the 
local  shock  Mach  number  (Eq.  28) 
according  to  the  Rankine-Hugoniot  condi¬ 
tion: 

p{o,y,t)/p]  =  1  +  [m2(y,t)-  1|2y/0'+  1) 

(30) 

This  produced  a  +  2  percent  pressure 
perturbation  near  the  wall. 

4.3  RESULTS. 

In  a  preliminary  calculation,  the  inflection 
point  was  placed  on  the  wall  (y  =  o).  After 
2000  cycles,  the  shear  layer  did  not  roll  up: 
instead,  it  merely  transported  the  initial 
shear  layer  profile  through  the  grid.  Hence, 
this  numerical  solution  satisfies  the 
Rayleigh-Tollmien  inflection  point  criterion 
which  states  that  shear  layers  without  an 
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inflection  point  are  inviscidly  stable  to 
small  perturbations. 

The  results  of  the  final  calculation  are 
depicted  in  Figure  10.  Frame  (a)  of  that 
figure  presents  the  overpressure  contours 
near  the  left  boundary  of  the  mesh  (0  <  x 
=$  200).  Acoustic  waves  radiate  from  the 
foot  of  the  shock,  similar  to  those  found  in 
the  experiment  (Figure  8).  Frames  (b)-(f) 
depict  the  rollup  of  the  shear  layer  into 
large  rotational  structures.  The  first  cusp 
in  the  material  interface  appears  between 
x  =  120  and  150;  by  x  =  220  the  structure 
has  rotated  through  one-half  revolution, 
and  a  complete  revolution  is  visible  at  x  = 
280.  There  is  some  evidence  of  incipient 
vortex  pairing  near  x  =  410,  470  and  650, 
but  the  Mode  I  response  seems  to 
predominate.  By  the  end  of  the  grid,  the 
visual  spreading  rate  is  about  8 '  vjS  =  0.01 5, 
which  agrees  with  the  optical  data  of  Figure 
8. 

Figure  1 1  presents  a  detailed  view  of  the 
flow  near  the  end  of  the  fine-grid  (800  < 
x  «  1000).  The  material  interface  plot 
shows  that  the  vortex  cores  have  rotated 
through  many  revolutions.  The  pressure 
contours  indicate  that  there  is  a  low 
pressure  region  centered  on  each  vortex, 
and  a  recompression  region  between 
vortices.  The  smoothness  of  the  pressure 
contours  indicates  that  the  wall  shear  layer 
is  free  of  shocks  and  basically  isentropic. 
even  though  the  flow  is  transonic.  The 
vorticity  contours  show  that  vorticity  is 
concentrated  near  the  wall,  but  accumu¬ 
lated  into  the  vortex  structures.  There  is 
also  a  baroclinic  generation  of  vorticity 
near  the  ends  of  the  vortex  structures,  due 
to  the  entrainment  of  density  gradient 
through  the  pressure  gradient  of  each 
vortex. 

The  flow  field  along  the  wall  (y  *  0.5)  is 
shown  in  Figure  12.  The  velocities  and 


densities  strongly  fluctuate,  even  beyond 
the  range  of  the  wall  and  free-stream 
values.  The  vortices  created  significant 
pressure  oscillations  with  peak  values  of  as 
much  as  -22  to  +8  percent  of  P2- 

The  flow  field  time  histories  were  moni¬ 
tored  at  stations  x  «  400,  600.  800  and  950. 
These  were  integrated  in  time  over  the  last 
1000  cycles  of  the  calculation  (about  20 
periods  at  x  =  950^.  The  resulting  mean- 
flow  velocity  and  density  profiles  are 
depicted  in  Figure  1 3  after  transformation 
back  to  laboratory-fixed  coordinates  (sub¬ 
script  o  denotes  free-stream  conditions 
above  the  boundary  layer,  which  corre¬ 
spond  to  state  2  behind  the  shock).  The 
results  have  been  scaled  by  means  of  the 
boundary  thickness  8bl  (where  u  =  0.99 

U»): 


7  BL  *  y/ d  BL  <31> 

with  8Bl  -  7.6 5,  9.15,  10.11.  11.92  and 
12.89  for  x  =  0,  400.  600.  800  and  950. 
respectively.  This  figure  shows  the  change 
or  evolution  of  the  mean-flow  velocity  from 
the  inflow  profile  (x=0)  which  contains  an 
inflection  point,  through  a  transition  region 
(x=400  to  600)  where  the  memory  of  the 
inflow  profile  is  fading,  and  finally  to  the 
asymptotic  profile  (x=800  and  950).  This 
scaling  seems  to  collapse  the  calculated 
profiles  quite  well  in  the  outer  region  (0.5 
<  Obl)-  Near  the  wall,  the  velocity  in  the 
first  cell  increases  with  increasing  x.  This 
is  caused  by  the  vortex  structures  which 
induce  an  enhanced  velocity  near  the  wall 
(even  though  the  no-slip  condition  is 
satisfied  in  the  initial  conditions  and  in  the 
boundary  condition  on  the  left  of  the  grid). 
Similar  comments  apply  to  the  mean 
density  profiles.  The  vortices  entrain  hot. 
low  density  fluid  near  the  wall,  which 
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Calculated  flow  field  of  the  development  of  a  wall  shear  layer  behind  a  shock 
(t  =  3075):  (a)  Overpressure  contours  showing  the  radiation  of  acoustic  waves 
from  the  foot  of  the  shock  (x=y=o):  (b)-(f)  material  interface  plots  demonstrating 
the  rollup  of  the  shear  layer. 
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"igure  11.  Expanded  view  of  the  flowfield  from  Figure  10  near  the  end  of  the  gut  •  a 
material  interface,  showing  the  rollup  of  the  layer:  (b)  density  contours.  (c> 
overpressure  contours  (solid  lines  correspond  to  positive  values):  (d)  vorticity 


contours. 


i.08 


£<-0.04 
|  -0-08. 
-0.12 
-0.16 
-0.20 
-0.24 

i 


0  200  400  600  800  1000 

X 


0  200  400  600  800  1000 

X 


Figure  12.  Flow  field  along  the  wall  of  the  shock  tube  wall-shear-layer  calculation  at  t  = 
3075:  (a)  density;  (b)  streamwise  velocity;  (c)  overpressure:  (d)  dynamic 
pressure. 


causes  the  density  in  the  first  cell  to 
decrease  with  increasing  x,  Clearly,  to 
continue  to  satisfy  the  no-slip  boundary 
condition  on  the  wall  (i.e..  to  accurately 
calculate  the  laminar  sublayer),  a  viscous 
flow  treatment  with  very  fine  zones  is 
needed  near  y=o.  Nevertheless,  the  profiles 
look  qualitatively  like  a  boundary  layer. 

Figure  14  presents  the  fluctuating  flow 
profiles  corresponding  to  Figure  13.  The 
streamwise  velocity  fluctuations,  u ' .  peak 
near  the  wall  and  at  t|bl  =  1  (no  doubt, 
above  and  below  each  vortex  structure). 
The  transverse  velocity  fluctuations,  v ' . 
reach  a  maximum  value  at  pbl  0.6 
(approximately  at  the  height  of  the  vortex 
centers).  Note  that  the  fluctuating  compo¬ 
nents  of  the  flow  do  not  approach  zero  at 
t)bl  =  T.  in  fact,  they  extend  out  3  or  4 
boundary  layer  thicknesses.  The  fluctuating 
shear  stress  u'v  seems  qualitatively 

,  .  .  .  rill 

correct,  (e  g.,  the  sign  is  opposite  to  —  ). 

r)V 

and  the  peak  value  increases  from  station 
1  to  2.  and  then  decays  at  stations  3  and 
4.  The  peak  values  of  the  fluctuating  flow 
components  are  listed  in  Table  1.  The  peak 
density  fluctuations  agree  very  well  with 
measurements  of  a  hypersonic  boundary 
layer  (Owen  et  al..  1975).  The  peak  shear 
stresses  are  somewhat  smaller  than  the 
low-speed  shear  layer  results. 

Next,  we  examine  whether  the  results  of 
this  wall  shear  layer  calculation  are 
consistent  with  known  properties  of  turbu¬ 
lent  boundary  layers.  According  to  exper¬ 
imental  measurements,  turbulent  boundary 
layers  may  be  divided  into  three  regions: 
(1)a  laminar  sublayer,  where  0  ^  riBL  < 
-0.02:  (2'  a  wall  region,  covering  -0.02  < 
qgL  -  -0.2:  and  (3)  a  wake  region 

extending  over  -0.2  <  pgL  <  I  Viscous 
effects  play  a  dominant  role  in  the  laminar 


sublayer.  In  the  wall  region,  the  mean-flow 
velocity  scales  with  the  wall  shear  velocity 
uT,  and  has  the  following  logarithmic 
dependence  on  y  (White.  1974): 

u/ur  =  —  ft)  ( vuT/v)  +  B  ( 32 ) 

K 

where 

k  =  0.40 
B  =  5.5 

Note  that  the  slope  of  the  profile  in  these 
coordinates  is  related  to  von  Karman  s 
constant,  k.  This  equation  may  be  extended 
to  cover  the  wake  region  by  adding  a 
so-called  “wake  function."  Coles  (1956) 
found  that  a  sin2  wake  function  gave  the 
best  correlation  of  a  variety  of  turbulent 
boundary  layer  data  (White.  1974);  hence. 
Eq.  (32)  becomes: 

it/iiT  =  —  fij  (yuT/ v)  +  B  +  A  sin;  (  ) 

(33) 

where 

'2.35  for  flatplate  (Schultz-Grunov . 
AS  .  1940) 

0.65  for  channel  flow  (Laufer.  1951) 

This  formulation  could  not  be  used  directly 
because,  for  example,  the  fluid  viscosity 
does  not  exist  in  the  present  inviscid 
calculation.  For  convenience,  then  Eq.  (33) 
was  rescaled  in  terms  of  the  free-stream 
velocity  and  t|bl-  and  the  viscosity  term 
was  absorbed  into  constant  B  .  yielding: 

Ti/U^  =  m  tr\  [,jbl)  +  B'+A'  sin :^'IBL) 


(34) 
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Figure  14.  Time-averaged  fluctuating  flow  profiles  from  the  shock  tube  wali-shear-laye 
calculation:  (a)  streamwise  velocity;  (b)  transverse  velocity;  (c)  shear  stress 
(d)  density;  (e)  pressure.  (Symbols  □.  O,  A  and  +  denote  stations  at  \  =  400 
600.  800  and  950.  respectively.). 
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Figure  14.  Time-averaged  fluctuating  flow  profiles  from  the  shock  tube  wall-shear-layer 
calculation:  (a)  streamwise  velocity;  (b)  transverse  velocity:  (c)  shear  stress: 
(d)  density:  (e)  pressure.  (Symbols  □.  O,  A  and  +  denote  stations  at  x  =  400. 
600.  800  and  950.  respectively.)  (Concluded). 
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where 

m  =  ur/UooK  (35) 

B'  =  |—  (t}(duT/v)  +  B\uT/Uao  (36) 

K 

A'  =  AUx/U  oe  (37) 

Note  that  because  of  the  rescaling,  the 
profile  parameters  (m,  B'  and  A' )  depend 
on  the  wall  shear  velocity  which  can  vary 
with  x.  Exclusive  of  the  viscous  sub-layer, 
Eq.  (34)  should  describe  the  mean-flow 
velocity  profile  over  essentially  the  entire 
height  of  the  boundary  layer  (~0.02  <  hbl 
<  ]). 

For  comparison  purposes,  the  calculated 
mean-flow  velocity  profiles  are  replotted  in 
Figure  15  using  semi-logarithmic  coordi¬ 
nates.  In  the  wall  region  (tibl  <  0.2),  the 
profiles  fan  out  with  a  definite  dependence 


on  distance  behind  the  shock.  The  calcu¬ 
lated  points  were  fit  in  this  region  with 
straight  lines  (the  solid  curves)  shown  in 
Figure  15,  which  have  the  logarithmic 
slopes  m  listed  in  Table  3. 

Near  the  wall,  it  was  assumed  that  the  total 
shear  stress  was  constant  and  equal  to  its 
peak  value,  hence  the  shear  velocity  uT  can 
be  evaluated  from  the  definition: 


=  I  \u'  v'|max/^ 


00 


(38) 


This,  in  turn,  can  be  used  to  evaluate  the 

i 

local  skin  friction  coefficient,  c  f ,  according 


to: 


C'=  21-^1  Jfl. 

/  l  PoO  U  00 


(39) 


Table  3.  Boundary  Layer  Parameters  for  the  Shock  Tube  Case. 


Wall  Values 

Profile  Parameters 

Cf 

X 

&BL 

-U  '  V '  /ui 
00 

uT/U 

1  00 

(io-3) 

m 

B' 

A' 

K 

0 

400 

7.65 

9.15 

0.0013 

0.0360 

3.5 

0.106 

0.3691 

0.6209 

0.338 

600 

10.11 

0.0016 

0.0400 

4.0 

0.0910 

0.3984 

0.5916 

0.439 

800 

11.92 

0.0012 

0.0346 

2.9 

0.0755 

0.4336 

0.5564 

0.456 

950 

12.89 

0.00059 

0.0243 

1.4 

0.0628 

0.4343 

0.5557 

0.387 

r  the  shock  tube  vvall-shear-layer  calculation,  demonstrating  the  Law  of  the  Wall 
he  wake  region  (dashed  curve.  Eq.  34)  of  the  profiles.  (Symbols  □,  O,  A  and  + 
)0.  800  and  950,  respectively.)  The  double-chain-dashed  curve  denotes  the  n  =  7 
at  plate.  The  triple-chain-dashed  curve  represents  a  boundary  layer  profile  with  an 
measured  by  Schubauer  and  Spangenberg  (1960). 


The  local  skin  friction  coefficients  eva¬ 
luated  from  this  calculation  are  listed  in 
Table  3:  they  vary  between  values  of  1.4 
to  4  x  10'3  which  are  in  excellent  agreement 
with  smooth-wall,  turbulent  boundary  layer 
measurements  (1.5  to  3  x  10‘3)  of  Martin 
(1958). 

Given  the  logarithmic  slope  and  the  shear 
velocity,  one  can  solve  Eq.  (35)  for  the  von 
Karman  constant: 

k  =  Ut/U^  m  (40) 

The  “von  Karman  constants”  so  evaluated 
from  the  calculated  velocity  profiles  are 
listed  in  Table  3.  They  vary  between  values 
of  0.34  to  0.46  (with  a  mean  value  of  0.405) 
—  which  is  in  excellent  agreement  with  the 
published  values  of  0.40  and  0.41.  In  other 
words,  the  calculated  mean-flow  velocity 
profiles  (specifically,  the  combination  of 
the  slopes  of  these  curves  and  their  shear 
velocities)  are  consistent  with  logarithmic 
profiles  associated  with  turbulent  boundary 
layers  in  the  wall  region. 

Figure  15  also  shows  that  in  the  outer 
region  of  the  boundary  layer,  the  calculated 
points  tend  to  converge  to  a  single  profile 
|i.e..  u/U*>  =  f(T)BL)l  which  is  basically 
independent  of  distance  behind  the  shock 
front.  The  dashed  curve  in  that  figure 
represents  Eq.  (34)  for  x  =»  950  (constants 
A’  and  B’  were  evaluated  at  t(bl  =  1  and 
0.1 ).  The  calculated  profiles  have  the  same 
shape  as  the  wake  curve,  with  the  maxi¬ 
mum  differences  being  less  than  8  percent. 
This  agreement  is  quite  good,  considering 
that  the  sin2  function  is  an  approximation 
to  a  variety  of  incompressible  boundary 
layer  data.  In  other  words,  the  calculated 
mean-flow  velocity  profiles  are  consistent 


with  the  wake  profiles  associated  with 
turbulent  boundary  layers. 

According  to  Schlichting  (1968).  the 
mean-flow  velocity  profile  associated  with 
turbulent  boundary  layers  can  be  approxi¬ 
mated  by  a  pow'er  function: 

u/Ux.  0.99(1,  BL)'I"  (41) 

where  t)bl  =  y/&BL-  The  value  n  =  7 
correlates  turbulent  boundary  layer  data  for 
flat  plates  with  zero  pressure  gradient.  This 
profile  is  depicted  in  Figure  15  as  a 
double-chain-dashed  curve.  It_gives  consid¬ 
erably  fuller  profiles  (e.g.,  u/U»  *  0.6) 
near  the  wall  than  the  present  calculation 
(u/U *  -  0.1  to  0.2). 

Previously  Mirels  (1956)  reasoned  that  the 
turbulent  boundary  layer  behind  a  shock 
wave  should  be  similar  to  the  turbulent 
boundary  layer  on  a  flat  plate,  and  thus  he 
assumed  that  n  *  7  power  function  profile 
applied  to  the  shock  tube  case.  The  present 
wall  shear  layer  calculation,  however, 
indicates  that  mean-flow  velocity  profile  for 
the  shock  tube  case  does  not  follow  this  n 
=  7  power  function  —  instead,  it  resembles 
the  boundary’  layer  profiles  measured  for 
cases  with  adverse  pressure  gradients.  The 
triple-chain-dashed  curve  in  Figure  15 
represents  such  data  of  Schubauer  and 
Spangenberg  (1960),  and  it  can  be  seen  that 
present  calculations  are  very  similar  to  their 
measured  profile.  Hence,  it  is  logical  to 
ask:  Which  is  the  correct  profile  for  a 
boundary  layer  behind  a  shock  wave?  In 
other  words,  are  the  profiles  of  the  present 
results  an  artifice  of  the  calculational 
method,  or  are  they  real?  Accurate, 
experimentally-measured  profiles  (e.g..  by 
means  of  iaser-doppler  velocimetry)  could 
answer  this  question:  unfortunately,  such 
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data  is  not  available.3  Instead,  we  shall 
address  this  question  indirectly  by  consider¬ 
ing  a  boundary  layer  case  that  has  a 
plethora  of  velocity-profile  data  —  namely, 
the  tripped  flat  plate  case.  This  is  described 
in  the  following  section. 

4.4  TRIPPED  FLAT  PLATE  CASE. 

Consider  a  two-dimensional,  laminar 
boundary  layer  developing  along  a  smooth 
flat  plate  in  the  absence  of  any  external 
pressure  gradient.  The  Mach  number  of  the 
free-stream  flow  above  the  boundary  layer 
is  taken  as  M»  =  0.2,  so  that  the  flow  is 
nearly  incompressible.  Then,  assume  that 
a  geometric  protuberance  is  located  on  the 
smooth  wall.  This  protuberance  causes  two 
physical  effects:  (1)  it  creates  an  inflec¬ 
tion-point  in  the  streamwise  velocity 
profile:  and  (2)  the  Strouhal  shedding 
frequency  of  its  wake  perturbs  the  wall 
layer,  thus  leading  to  a  “turbulent” 
boundary  layer  flow. 

The  development  of  the  fluctuating  flow 
downstream  of  this  disturbance  was 
simulated  numerically  by  the  high-order 
Godunov  code.  The  computational  grid 
used  was  identical  to  that  of  the  shock  tube 
case.  The  flow  field  was  initialized  with  a 
uniform  density  (p  =  1)  and  pressure  (p  = 
17.86)  field.  The  wall  shear  layer  was 
modeled  by  a  Tanh(y)  streamwise  velocity 
profile: 

»(.v.  y.  o)/U<j>  =Av)  (42) 

=  0.5 1 1  +  Tunh(0.5y-  1.5)Jv(t,>\  o)  =  0 

Note  that  by  construction,  u  =  0  on  the  wall 
thereby  satisfying  the  no-slip  condition, 


while  in  the  free  stream  u  =  Uoo-  The 
pressure  was  chosen  so  that  the  free  stream 
Mach  number  was  Mx  =  01  The  inflection 
point  and  the  initial  tracer  particles  were 
located  at  y  =  3. 

The  left  boundary  of  the  grid  was  driven 
with  the  same  flow  field  with  a  sinusoidal 
perturbation  on  the  streamwise  velocity: 

u  (o,y,t)  /U  * 

10 

=  /(y)  (1+00')  X  €n  sin  u)„t\  (43) 
n  =  1 

<P  (y)  =  1 1  -/  (y)l  (44) 

with  € |  =0.10.  Note  that  by  means  of  Eq. 
(44),  the  perturbations  reach  a  maximum 
value  ( *  0.025  U  »)  at  about  half  the  height 
of  the  shear  layer  and  are  zero  outside  this 
region. 

The  calculation  was  run  104  cycles  to 
eliminate  transients  from  the  start-up 
process  and  to  build  up  data  for  analyzing 
the  mean-flow  profiles.  The  flow  field  at  the 
end  of  the  calculation  (t  =  361 7)  is  depicted 
in  Figure  16.  The  tripped  wall  shear  layer 
rolled  up  into  large  rotational  structures 
similar  to  the  shock  tube  calculation,  except 
in  this  case  the  structures  rotated  in  a 
clockwise  direction.  The  material  interface 
plots  (Fig.  16  a,b,c)  show  that  the  first  cusp 
appears  at  x  s  40,  and  the  vortex  structure 
has  rotated  through  one-half  a  revolution 
at  x  *  55.  By  x  *  90  the  vortices  are 
well-formed  and  have  rotated  through  one 
revolution.  Two  cores  are  seen  in  the 
vortices  at  x  =  210,  320  and  430  (i.e..  every 
other  vortex),  thus  indicating  a  Mode  II 
pairing,  but  multiple  merging  did  not  seem 


3.  Martin  (1958)  did  report  velocity  profiles  for  a  boundary  layer  behind  a  shock.  However,  these 
velocities  were  not  measured  directly  but  were  inferred  from  interferometric  measurements  of  density 
profiles  by  use  of  the  Crocco  relation.  This  relation  is  valid  only  along  streamlines  and  cannot  be 
applied  in  turbulent  boundary  layers  because  of  the  large-scale  rotational  structures. 
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Figure  16.  Flow  field  of  the  tripped  flat  plate  calculation  (t  =  3617):  (a),  (b).  (c)  are  material 
interface  plots  showing  the  rollup  of  the  layer;  (d),  (e)  and  (f)  represent  contour 
plots  of  density,  o'  erpressure  and  vorticity,  respectively. 


42 


to  occur  even  though  the  first  ten  pertur¬ 
bation  frequencies  were  used.  Also  shown 
in  Figure  16d,e.f  are  the  density,  overpres¬ 
sure  and  vorticity  plots  corresponding  to 
Figure  16c.  The  maximum  pressure  and 
density  variations  are  less  than  one-half 
percent  for  this  low  speed  flow. 

Clearly,  this  calculation  captures  only  the 
largest-scale  rotational  structures  that  are 
essentially  independent  of  Reynolds  num¬ 
ber.  Smail-scale  vortices,  which  should  be 
embedded  in  the  large  rotational  structures, 
and  which  are.  no  doubt,  Reynolds  number 
dependent  (e.g.,  as  found  in  the  Brown  and 
Roshko  free  shear  experiments),  are  not 
resolved  on  this  grid. 

The  flow  field  was  monitored  at  stations  x 
=  400.  600,  800,  950.  This  data  was 
integrated  over  the  last  4000  cycles  to 
establish  the  time-averaged  profiles.  The 
resulting  mean-flow,  streamwise  velocity 
profiles  are  presented  in  Figure  17.  The 
profiles  are  qualitatively  similar  to  the 
previous  calculation  in  that  they  converge 
to  a  single  profile  in  the  wake  region  and 
exhibit  a  logarithmic  behavior  in  the  wall 
region.  Now,  however,  the  magnitudes  of 
the  velocity  are  very  similar  to  those 
obtained  from  the  n  =  7  power  function 
profile  depicted  as  the  chain-dashed  curve. 
No  doubt,  the  flow  at  station  x  =  400  is 
transitional. 

The  calculated  profiles  were  fitted  with  the 
Coles  boundary  layer  function  (Eq.  34), 
and  the  profile  parameters  are  listed  in 
Table  4.  The  von  Karman  constants 
inferred  from  these  fits  (0.37  <  k  <  0.50) 
are  quite  reasonable.  Comparing  with 
available  experimental  data  (Table  4),  we 
find  that  local  skin  friction  coefficients  are 
too  small  by  a  factor  of  6,  and  the 
associated  shear  velocities  are  too  small  by 
a  factor  of  2.5.  Wieghardt  and  Tillman 
(1944)  have  performed  some  of  the  most 


detailed  and  extensive  measurements  of 
mean  velocity  profiles  for  a  tripped, 
turbulent  boundary  layer  on  a  flat  plate: 
this  set  of  data  is  compiled  in  Coles  and 
Hirst  (1969).  Two  of  these  stations,  x  = 
79cm  (profile  1400-10)  and  x  =  499cm 
(profile  1400-26),  which  cover  the  x/8 
range  of  the  numerical  simulation,  are 
shown  in  Figure  17.  The  dashed  curves  on 
that  figure  represent  the  Coles  bound- 
ary-layer-function  fits  to  that  data,  e.g.. 

tT  /«4  =  0.0779  In  nBL 

+  0.849  +  0.141  sin^  (y  i]qi)  (45) 

for  station  x  =  79cm.  This  data  may  be  used 
to  evaluate  the  accuracy  of  the  calculation. 

Comparing  the  solid  curve  of  the  calcula¬ 
tion  (x  =  600;  x/8  =  57)  with  the  lower 
dashed  curve  (Eq.  45)  representing  the 
Wieghardt  and  Tillmann  data  (x  =  79cm: 
x/8  *  56),  one  sees  that  the  two  curves  are 
identical— except  for  the  first  calculated 
point.  Similar  conclusions  apply  to  the  two 
downstream  stations;  outside  the  first  row 
of  celis.  the  calculated  velocities  agree  very 
well  with  the  data  (e.g..  less  than  4  percent 
difference).  In  the  first  row  of  cells, 
however,  the  calculated  velocities  are  7  to 
14  percent  too  large,  probably  because 
viscous  losses  near  the  wall  are  not 
included  in  the  calculation  (e.g.,  the  viscous 
sublayer  extends  to  approximately  20 
percent  of  the  height  of  the  first  row  of 
computational  cells,  and  large  momentum 
losses  in  this  sublayer  due  to  viscosity  could 
decrease  the  cell-averaged  momenta  in 
these  cells  by  as  much  as  20  percent). 

The  root-mean-squared  flow  fluctuations 
were  also  evaluated  for  this  case.  They  are 
compared  with  the  hot-wire  anemometry 
measurements  of  Klebanoff  (1955)  in 
Figure  18  (dashed  curves),  and  the  peak 
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Figure  17.  Mean-flow  velocity  profiles  of  the  tripped  flat  plate  case.  Symbols  □.  O.  A.  and  O  denote  calculation  stations 
at  x  =  400.  600,  800  and  950,  respectively.  The  chain-dashed  curve  represents  the  n  =  7  power  function  profile 
(Eq.  41).  The  dashed  curves  denote  the  Coles  function  fits  of  the  tripped  flat  plate  data  of  Wieghardt  and  Tillmann 
(1944):  symbols  •  and  ♦  denote  data  stations  at  x  =  79  cm  (profile  1400-10)  and  x  =  499  cm  (profile  1400-26), 
respectively. 
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Figure  18.  Time-average  fluctuating  flow  profiles  from  the  tripped  flat  plate  calculation: 

(a)  streamwise  velocity;  (b)  transverse  velocity;  (c)  shear  stress:  (Symbols  □. 
O.  A.  and  +  denote  stations  at  x  =  400,  600,  800  and  950.respectively>.  Dashed 
curves  represent  the  data  of  Klebanoff  (1955). 
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Table  4.  Boundary  Layer  Parameters  for  the  Tripped  Flat  Plate  Case. 


1  -  — 

Stations 

Wall  Values 

Profile  Parameters 

-u'v  7Ui 
00 

IVU. 

Cf' 

X 

8 

x/8 

(10-3) 

(10°) 

m 

B' 

A' 

K 

Calculation 

0 

7.65 

0 

- 

-- 

-- 

-- 

- 

-- 

- 

400 

9.5 

42 

0.531 

0.0230 

1.06 

0.0619 

0.763 

0.227 

0.37 

600 

10.5 

57 

0.259 

0.0161 

0.518 

0.0428 

0.789 

0.201 

0.38 

800 

11.5 

70 

0.226 

0.0150 

0.452 

0.0319 

0.797 

0.193 

0.47 

950 

12.5 

76 

0.186 

0.0136 

0.372 

0.0273 

0.806 

0.184 

0.50 

Wieghardt  &  Tillmann(1944) 

79cm 

1 ,41cm 

56 

1.59a 

0.0398 

3.17 

0.0779 

0.849 

0.141 

0.51 

499cm 

6.45cm 

77 

1.22“ 

0.0349 

2.43 

0.0783 

0.880 

0.110 

0.45 

Klebanoff  (1955) 

1.38“ 

0.0371 

2.75 

- 

- 

- 

- 

Schlichting  (1968) 
(Eq.  2l.12:Rex  = 

1.9x107) 

1.04“ 

0.0322 

2.07 

- 

- 

- 

- 

■ 

aEvaluated  from  Eq.  38 


values  are  listed  in  Table  1.  From  these 
comparisons  it  is  obvious  that  th  r.m.s. 
flow  fluctuations  are  predicted  quite  poorly 
by  the  present  calculation,  even  though  the 
mean  profiles  are  in  good  agreement  with 
the  data.  Presumably  this  is  because 
fine-scale  vortices,  which  should  be  em¬ 
bedded  within  the  large  rotational  struc¬ 
tures.  are  not  resolved  in  this  calculation 
and  because  three-dimensional  effects  are 
not  included. 

4.5  SUMMARY. 

We  conclude  that  this  tripped  flat  plate 
calculation  qualitatively  simulates  some  of 


the  features  found  in  the  wake  region  of 
turbulent  boundary  layers  on  flat  plates. 
The  calculated  mean-flow  velocity  profiles 
agree  with  the  data  of  Wieghardt  and 
Tillmann,  except  for  the  first  row  of  cells 
on  the  wall.  This  agreement  with  data 
implies  that  our  calculated  profiles  are 
neither  fortuitous  nor  an  artifice  of  the 
numerical  simulation,  and  that  the  calcula¬ 
tion  qualitatively  represents  the  largest 
fluid-dynamic  features  resolvable  on  the 
computational  grid. 

The  above  comparisons  offer  aposteriori 
evidence  that  the  velocity  profiles  of  the 
shock  tube  calculation  are  also  qualitatively 
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correct.  The  large-scale  features  of  the  flow 
field  (e.g.,  the  formation  of  discrete  density 
and  vortex  structures,  the  generation  of 
acoustic  waves  at  the  foot  of  the  shock 
front,  and  the  visual  spreading  rate)  agree 
with  the  available  shock  tube  data.  The 
mean-flow  velocity  profiles  exhibit  a  log(y) 
velocity  distribution  near  the  wall,  while  in 
the  outer  region  the  profiles  collapse  to  a 
single  curve  that  behaves  like  the  wake 
function  of  a  turbulent  boundary  layer. 
Hence,  it  appears  that  the  turbulent 
boundary  layer  profile  behind  a  shock  (e.g.. 
Fig.  15)  is  distinctly  different  from  the  flat 
plate  case  (e.g..  Fig.  17).  Perhaps  this 
difference  is  related  to  the  slope  of  the 
mean  velocity  profile  (fUT/fly ).  which  is 
negative  for  the  shock  tube  case  (in 
shock-fixed  coordinates)  and  positive  for 
the  flat  plate  case. 

Referring  back  to  Figures  1  la  or  16c.  it  can 
be  seen  that  the  logarithmic  velocity  region 
near  the  wall  is  simply  the  time-averaged 
velocity  field  induced  below  each  vortex, 
while  the  wake  profile  corresponds  to  the 
time-averaged  velocity  field  within  and 


above  each  vortex  structure.  In  other 
words,  the  fluctuating  flow  associated  with 
turbulent  boundary  layers  is  merely  the 
evolution  of  a  perturbed  but  unstable  wall 
shear  layer.  This  concept  is  not  new:  it  has 
been  expressed  by  Coles  (1969).  and 
others. 

These  wall-shear-layer  calculations  capture 
only  the  largest-scale.  inviscid  rotational 
structures.  Small-scale  vortices  (which  may 
depend  on  the  Reynolds  number)  will  be 
shed  from  the  wall  and  entrained  into  the 
large  rotational  structures,  but  a  much  finer 
computational  grid  would  be  required  to 
resolve  such  effects.  Viscous  forces  must 
be  included  to  correctly  model  the  laminar 
sublayer.  Eventually  the  boundary  layer 
will  develop  three-dimensional  vortex 
structures,  and  accurate  numerical  simula¬ 
tions  must  also  include  these  effects. 
Nevertheless,  the  above  agreement  with 
data  indicates  that  such  features  as  the 
mean-flow  velocity  profiles  are  dominated 
by  two-dimensional  inviscid  effects  over 
most  of  the  boundary  laver  (e.e..  0.1  <  riBL 
<  1). 
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SECTION  5 


WALL  JET  CALCULATION 


The  final  case  considered  is  a  numerical 
simulation  of  a  low-speed  wall  jet  propagat¬ 
ing  into  a  quiescent  fluid,  similar  to  the 
water  tunnel  experiments  of  Bajura  and 
Catalano  (1975).  The  calculation  dramati¬ 
cally  demonstrates  the  strong  coupling 
between  the  dynamics  of  the  free  shear 
layer  on  top  of  the  jet  and  the  wall  boundary 
layer  below  the  jet. 

5.1  FORMULATION. 

The  calculation  was  performed  on  a  grid 
with  a  fine-zoned  region  of  200  x-cells  by 
80  y-cells.  (Ax=Ay=1).  as  shown  schemati¬ 
cally  in  Figure  19.  Expanding  cells  were 
used  above  and  to  the  right  of  the 
fine-zoned  region,  resulting  in  a  total 
domain  of  8500  by  1500.  An  inviscid  (i.e.. 
slip  wall)  boundary  condition  was  em¬ 
ployed  along  the  wall,  while  an  outflow 
boundary  condition  was  used  at  the  right 
edge  of  the  grid. 

The  flow  field  was  initialized  with  a  wall 
jet  velocity  profile  consisting  of  a  Tanh(v) 
profile  for  the  free  shear  layer  portion  and 


a  Polhausen  profile  representing  (he  wall 
boundary  layer: 

u  {x,y,o) 

0.5  Uj  fl  -  Tanh(y~  1 7) |  for  y  >  5 
=  '  2  Uj  £  (1  - 1  2  +  0.5  £  3  |  for  v  <  5  ’ 

*  d 

v(x,y,o)  =  0  1 46) 

where 

Uy  =  2 

I »  y/to 

The  initial  shear  thicknesses  were  8„  =  5 
and  6  on  the  wall  and  free  shear  layers, 
respectively.  The  center  line  of  (he  free 
shear  layer  was  at  y  *  17  (=  H).  The  jet 
was  initialized  with  uniform  thermodynam¬ 
ic  properties: 

p  ==  1 ;  p  =  71 .4:  e-  179. 

The  pressure  was  selected  so  that  the  peak 
value  of  the  Mach  number  of  the  jet  was 
Mj  =  0.2. 


Grid 
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The  left  boundary  of  the  grid  was  driven 
with  the  same  wail  jet  profile,  using  a 
sinusoidal  perturbation  on  the  streamwise 
velocity  component  (Eq.  15).  The  first  ten 
modes  were  employed,  with  a  maximum 
perturbation  amplitude  of  £|  =  0.01. 
Material  lines,  initially  located  at  the 
inflection  point  of  the  free  shear  layer  (/ 
=  17)  and  in  the  boundary  layer  (y  =  2), 
were  tracked  to  follow  the  interactions  of 
the  two  shear  layers. 

5.2  RESULTS. 

The  calculation  was  run  8000  cycles.  Figure 
20  presents  a  series  of  frames  of  the 
calculated  material  interface  shapes  cover¬ 
ing  one  complete  period  during  the 
calculation.  The  merging  patterns  are  quite 
different  from  those  found  in  the  above-de¬ 
scribed  shear  layer  calculations  (due  to 
differences  in  \.  r  and  \p).  In  general,  the 
free  shear  layer  becomes  unstable  first4, 
and  forms  positive  vortex  structures 
(labeled  sequentially  with  cardinal  num¬ 
bers).  Induced  velocities  from  these  struc¬ 
tures  cause  the  wall  shear  layer  to  become 
unstable,  and  form  negative  vortex  struc¬ 
tures  (labeled  sequentially  with  capital 
letters),  that  explode  off  the  wall.  This 
behavior  is  qualitatively  similar  to  that 
found  in  the  wall  jet  experiments  of  Bajura 
and  Catalano  (1975).  Large-scale  struc¬ 
tures  from  the  two  shear  layers  interact  and 
pair,  thus  creating  a  von  Karman-type 
vortex  street  of  alternating  positive  and 
negative  vortex  structures. 

The  following  is  a  detailed  description  of 
this  process.  Frames  (a)  through  (g)  of 
Figure  20  show  the  formation  of  structure 
1 -2-3-4.  Vortices  1  and  2  pair  and  merge. 


and  induce  the  rollup  of  vortex  A  on  the 
wall  (frames  a  and  b).  Structure  1-2 
entrains  vortex  3  from  below  (frames  c  and 
d).  Vortex  4  triggers  the  creation  of  vortex 
B  on  the  wall  (frames  d  and  e).  Then 
structure  1-2-3  entrains  vortex  4  from 
below,  which  causes  the  merger  of  wall 
vortices  A-B  (frames  e  and  f).  Finally,  the 
large-scale  structure  1 -2-3-4  interacts  with 
the  structure  A-B.  causing  it  to  explode  off 
the  wall  (frame  g).  These  structures 
(1 -2-3-4  and  A-B)  then  interact  to  form  the 
first  pair  of  a  von  Karman-type  vortex  street 
which  passes  off  the  grid.  The  creation  of 
structure  1-2-3-4  is  here  labeled  a  Sequen¬ 
tial  Merger  of  four  vortices. 

Meanwhile,  vortices  5  and  6  have  paired, 
inducing  the  rollup  of  vortex  C  on  the  wall 
(frames  d  and  e);  also  vortices  7  and  8  hav  e 
paired,  inducing  the  rollup  of  vortex  D  on 
the  wall  (frames  f  and  g).  Then  the  pair  5-6 
merges  with  pair  7-8.  forcing  wall  vortices 
C-D  to  merge  (frames  h  and  i).  Finally,  the 
large-scale  structure  5- 6-7-8  interacts  with 
structure  C-D,  causing  it  to  explode  off  the 
wall  (frames  i  and  j).  These  structures 
(5-6-7-8  and  C-D)  then  interact  to  form  the 
second  pair  of  a  von  Karman-type  vortex 
street  (frames  k  and  I)  which  passes  off  the 
grid.  The  creation  of  the  large-scale 
structure  5-6-7-8  is  here  labeled  a  Paired 
Merger  of  four  vortices. 

This  Paired  Merger  process  is  repeated 
creating  large-scale  structures  9-10-11-12 
and  E-F  which  form  the  third  pair  of  a  von 
Karman-type  vortex  street  (frames  o  and  p) 
that  passes  of  the  grid. 

This  is  followed  by  a  Sequential  Merger  of 
vortices  13,14,15  and  16  (frames  I  through 
p)  and  the  creation  of  the  wall  structure 


4.  *  Preliminary  calculations  used  a  Tanli(y)  profile  with  an  inflection  point  off  the  wall  in  place  nl  a 

Polhausen  profile.  These  results  also  showed  that  the  free  shear  layer  became  unstable  first.  Hence,  the 
results  were  insensitive  to  the  profile  used  for  the  wall  shear  layer. 
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Figure  20.  Numerical  simulation  of  a  wall  jet  (\  «  -1,  r  =  oo,  =  0).  The  figures  depict 
the  evolution  of  the  material  interfaces  over  one  period  of  the  flow.  Vortices 
forming  on  the  free  shear  layer  are  labeled  sequentially  with  cardinal  numbers, 
while  those  forming  on  the  wall  layer  are  denoted  by  capital  letters. 
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Figure  20.  Numerical  simulation  of  a  wall  jet  (K  =  -1,  r  oo,  \p  =  0).  The  figures  depict 
the  evolution  of  the  material  interfaces  over  one  period  of  the  flow.  Vortices 
forming  on  the  free  shear  layer  are  labeled  sequentially  with  cardinal  numbers, 
while  those  forming  on  the  wall  layer  are  denoted  by  capital  letters  (Continued). 
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Figure  20.  Numerical  simulation  of  a  wall  jet  (K  =  -1,  r  =  oo,  \p  =  0).  The  figures  depict 
the  evolution  of  the  material  interfaces  over  one  period  of  the  flow .  V  ortices 
forming  on  the  free  shear  layer  are  labeled  sequentially  with  cardinal  numbers, 
while  those  forming  on  the  wall  layer  are  denoted  by  capital  letters  (Continued). 


52 


0  50  100  150  200 

X 


Figure  20.  Numerical  simulation  of  a  wall  jet  (V  *  -1,  r  *  oo,  \p  =  0).  The  figures  depict 
the  evolution  of  the  material  interfaces  over  one  period  of  the  flow.  Vortices 
forming  on  the  free  shear  layer  are  labeled  sequentially  with  cardinal  numbers, 
while  those  forming  on  the  wall  layer  are  denoted  by  capital  letters  (Concluded). 
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G-H.  These  structures  interact  to  form  the 
fourth  pair  of  a  von  Karman  vortex  street 
(frames  s  and  t)  that  passes  off  the  grid. 

Meanwhile,  vortices  17  and  18  have  paired 
(frames  p  and  q)  creating  wall  vortex  1; 
vortices  19  and  20  have  paired  (frames  r 
and  s)  creating  wall  vortex  J;  and  vortices 
1  and  2  have  paired  (frame  t)  creating  wall 
vortex  A.  Note  that  the  interface  shape  in 
frame  (t)  is  essentially  identical  to  that  of 
frame  (a).  Hence,  the  whole  process  is  then 
repeated.  Also  note  that  the  vortex  pair 
17-18  and  the  pair  19-20  pass  off  the  grid 
without  merging  (frames  a  through  f). 

In  summary,  the  vortex  merging  on  the  free 
shear  layer  of  the  wall  jet  follows  this 
sequence: 

•  1  Sequential  Merger  (1 -2-3-4) 

•  1  Paired  Merger  (5 -6- 7-8) 

•  1  Paired  Merger  (9-10-11-12) 

•  1  Sequential  Merger  (13-14-15-16) 

•  1  Pairing  (17-18) 

•  1  Pairing  (19-20) 

and  then  the  sequence  repeats.  The 
Sequential  Mergers  correspond  to  a  Mode  IV 
response,  while  the  Paired  Mergers  corre¬ 
spond  to  two  Mode  11  responses  followed 
by  a  Mode  IV  response.  Evidently  the 
higher  modes  did  not  directly  couple  into 
the  dynamics  for  these  flow  lengths. 

Figure  21  depicts  the  flow  field  contours 
corresponding  to  the  last  frame  of  Figure 
20.  The  vorticity  on  the  free  shear  layer 
rapidly  accumulates  into  the  large  struc¬ 
tures  as  a  result  of  the  first  pairing  of 
vortices,  and  thus,  the  braid  regions  are 
essentially  devoid  of  vorticity.  A  low 
pressure  region  is  formed  in  the  center  of 


each  vortex  structure  due  to  the  rotational 
flow,  while  a  recompression  is  created  in 
the  braid  region  between  structures.  The 
vorticity  in  the  wall  shear  layer  accumulates 
but  it  remains  near  the  wall  until  the  wall 
vortex  structures  are  entrained  into  the  von 
Karman  vortex  street. 

The  flow  field  time-histories  were  moni¬ 
tored  at  stations  x  =  85  and  170  (corre¬ 
sponding  to  5H  and  10H).  These  were 
integrated  in  time  over  the  last  4000  cycles 
of  the  calculation  (with  about  8  large-scale 
structures  passing  station  170).  The  result¬ 
ing  mean  velocity  profiles  are  presented  in 
Figure  21.  The  results  have  been  scaled  by 
the  jet  thickness  8j  (where  u  =  0.5  um;,x): 

rjj^y/dj  (47) 

where  8j  =  17,  21.7  and  31.2  for  \  =  0.  85 
and  170,  respectively.  Figure_22  demon¬ 
strates  that  the  mean  velocity  u  profiles  at 
x  *  85  and  170  are  considerabb  different 
than  the  inflow  profile  (shown  as  a  solid 
curve).  The  peak  value  of  u  decays  with 
distance,  indicating  that  the  streamwise 
momentum  is  spread  lateral ly  due  to 
large-scale  mixing.  The  shear  layer  portion 
of  the  profiles  (0.5  <  qj)  at  the  two  stations 
collapses  quite  well  with  this  scaling.  The 
wall  layer  portion  of  the  u  profiles  (r)j  <  0.5) 
is  significantly  modified  by  the  vortex 
structures  of  the  free  shear  layer.  By  \  = 
170,  a  separated  flow  profile  has  devel¬ 
oped.  The  large-scale  structures  induce  a 
mean  transverse  velocity  v  across  the  entire 
width  of  the  grid. 

Figure  23  presents  the  fluctuating  flow 
profiles  (corresponding  to  Figure  22).  while 
the  peak  values  are  listed  in  Table  1.  The 
streamwise  velocity  fluctuations  peak  near 
the  wall  (t|j  *»  0.2)  and  near  the  free  shear 
layer  (rij  *  0.75),  similar  to  the  results  of 
Bajura  and  Catalano  ( 1975).  The  calculated 
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Figure  21.  Flow  field  of  the  wall  jet  calculation  at  t  =  1722:  (a)  material  interfaces;  (b) 
vorticity  contours  (solid  lines  correspond  to  positive  values,  while  dashed  lines 
denote  negative  values):  (c)  overpressure  contours  (solid  lines  correspond  to 
positive  values). 
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Figure  23.  Time-averaged  fl'.^uating  flow  profiles  from  the  wall  jet  calculation:  (a)  stream- 
wise  velocity:  (b)  transverse  velocity:  (c)  shear  stress.  (Symbols  □  and  O  denote 
stations  at  x  =  85  and  170.  respectively.). 
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peaks,  however,  are  an  order  of  magnitude 
larger  than  the  data,  presumably  due  to  the 
much  larger  shear  used  in  the  calculation. 
(Their  experiments  had  an  exit  Reynolds 
number  of  only  445.)  The  transverse 
velocity  fluctuations  also  have  two  peaks, 
one  at  t|j  =  0.5  and  the  other  at  t|j  =  1,  at 
the  last  station.  The  fluctuating  shear 
stresses  are  positive  and  have  a  single  peak 
at  r)j  =  0.75,  indicating  that  the  structures 
on  the  shear  layer  are  controlling  the 
mixing  process. 

5.3  SUMMARY. 


The  dynamics  of  the  wall  jet  was 
considerably  different  than  the  previously 
described  cases.  The  free  shear  layer  on  top 
of  the  jet  reacted  first  by  forming  discrete 
vortex  structures  which  paired  and  merged. 
The  merging  patterns  consisted  of  both 
paired  and  sequential  merging  of  four 
vortices,  corresponding  to  both  Mode 
ll-type  and  Mode  IV-type  response  —  even 
though  the  first  ten  subharmonics  were 
used  to  perturb  the  jet. 


The  merging  was  faster  than  in  the 
spreading  shear  layer  calculation  because 
of  the  larger  absolute  shear  and  it  was 
smoother  and  more  periodic  because  the 
uniform  density  assumption  of  the  wall  jet 
calculation  virtually  eliminated  am  baro- 
clinic  generation  of  vorticity.  Vortex 
pairing  on  the  free  shear  layer  triggered 
instabilities  in  the  wall  shear  layer  which 
rolled  up  into  a  vortex  directly  underneath 
each  large  structure  on  the  free  shear  layer. 
Merging  of  these  latter  structures  caused 
the  wall  vortices  to  merge.  In  the  final  stage 
of  interactions,  the  large-scale  structures 
from  the  two  shear  layers  paired,  forming 
a  von  Karman-type  vortex  street.  In  this 
way,  fluid  from  the  wall  layer  of  the  jet  was 
mixed  across  the  free  shear  layer,  and  vice 
versa.  In  effect,  the  rotational  structures  on 
the  free  shear  layer  caused  the  wall  layer 
to  periodically  explode  off  the  wall,  thus 
creating  a  separated  flow  profile  These 
results  are  qualitatively  similar  to  the  wall 
jet  experiments  of  Bajura  and  Catalano 
(1975).  Peak  velocity  fluctuations,  howev¬ 
er.  are  larger  in  this  calculation  because  of 
the  faster  merging  and  larger  absolute 
shear. 
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SECTION  6 


CONCLUSIONS  AND  GENERALIZATIONS 


The  present  inviscid  calculations  seem  to 
duplicate  quite  well  the  major  flow  features 
that  have  been  observed  in  experiments. 
These  include  the  formation  and  growth  of 
large  rotational  structures,  the  visual 
spreading  rates,  and  the  mean-flow  pro¬ 
files.  Presumably,  these  features  are  domi¬ 
nated  by  two-dimensional  effects.  The 
fluctuating  components  of  the  flow  agree 
qualitatively  with  the  available  experimen¬ 
tal  data,  hpwever.  peak  fluctuations  such 
as  u'  or  v’  can  be  1.5  to  2  times  larger 
than  the  data  because  of  the  two-dimen¬ 
sional  flow  approximation. 

It  is  now  clear  that  the  fluctuating,  time- 
dependent  flow  variations  associated  with 
these  “steady"  shear  layer  problems  are 
caused  by  the  dynamic  evolution  of 
unstable  shear  layers.  By  taking  into 
account  one  more  degree  of  freedom  (time) 
in  these  “steady"  shear  layer  calculations, 
one  can  determine  not  only  the  mean  flow 
but  also  the  fluctuating  flow  field  without 
any  turbulence  modeling.  Since  the  present 
calculations  are  in  good  agreement  with  the 
available  experimental  data,  we  conclude 
that  the  evolution  of  these  unstable  shear 
layers  is  dominated  by  inviscid  flow  effects. 
The  calculations  also  point  out  the  im¬ 
portance  of  accurately  evaluating  the 
convective  derivatives.  Viscous  forces  and 
molecular  diffusion  are  too  small  (relative 
to  convective  effects)  to  significantly  affect 
the  dynamics  of  •  such  large-Reynolds- 
number  flows.  Molecular  effects  are  then 
relegated  to  the  relatively  minor  role  of 
de-singularizing .  the  inviscid  problem  by- 
spreading  the  shear  over  a  small  but  finite 
volume  of  fluid,  and  they  create  additional 
vorticity  at  wall  boundaries. 


In  fact,  molecular  transport  effects  such  as 
fluid  diffusivity  and  viscosity  are  inherently 
diffusional  effects.  They  correspond  to 
fluctuations  on  a  molecular  scale,  but  on 
a  macroscopic  scale  they  can  only  produce 
smoothing  —  that  is,  they  can  never  by 
themselves  create  the  self-organized  behav¬ 
ior  such  as  the  large-scale  rotational 
structures  presented  here.  Consequently, 
the  macroscopic  fluctuations  of  the  flow  are 
caused  exclusively  by  inviscid.  convective 
flow  effects  —  that  is,  by  the  nonlinear 
dynamics  of  unstable  shear  layers. 

The  calculations  presented  here  show  that 
the  fluid-dynamic  response  of  three  funda¬ 
mentally  different  shear  layers  is  quite 
similar  in  that  they  roll  up  into  large-scale 
rotational  structures  which  grow  by  interac¬ 
tions  and  merging.  There  are  five  features 
which  these  shear  layers  have  in  common, 
namely: 

1.  The  flows  contain  an  initial  tangen¬ 
tial  velocity  change  (i.e..  a  shear 
layer). 

2.  The  initial  shear  layer  thickness  is 
non-zero  (i.e..  the  problem  is 
desingularized). 

3.  The  velocity  profiles  are  unstable  to 
infinitesimal  perturbations. 

4.  Perturbations  exist  in  the  flow 
(either  from  external  sources,  or 
from  internal  sources  such  as  the 
growth  of  molecular  fluctuations). 

5.  The  Reynolds  number  of  the  flows 
are  large  (i.e..  the  convective  flow 
velocities  are  orders  of  magnitude 
larger  than  any  molecular  diffusion 
velocities). 


These  five  common  features  seem  to  be 
responsible  for  the  development  of  the 
organized  rotational  structures  in  these 
(and  probably  other)  shear  layers. 

Successful  numerical  simulations  of  such 
fluctuating  flows  have  numerical  require¬ 
ments  which  parallel  the  above  features, 
namely: 

1.  They  contain  the  correct  amount  of 
circulation,  by  means  of  initial  and 
boundary  conditions. 

2.  The  initial  shear  layer  is  resolved  on 
the  computational  mesh. 

3.  The  initial  velocity  profile  contains 
an  inflection  point,  hence  the  flow 
is  unstable. 

4.  The  shear  layers  are  perturbed  with 
a  spectrum  of  frequencies. 

5.  The  numerical  Reynolds  number  of 
the  finite  difference  solution  is  large 
so  that  convective  flow  effects  are 
accurately  evaluated  (i.e..  the  algo¬ 
rithm  is  nondiffusive  enough  to 
allow  the  instabilities  to  grow,  and 
not  be  artificially  damped  by  numer¬ 
ical  diffusion). 

Accurate  numerical  simulations  have  two 
additional  requirements: 

6.  The  calculations  must  be  at  least 
two-dimensional  and  time  depen¬ 
dent.  even  for  problems  which  are 
steady  in  the  mean  or  time-averaged 
sense. 

7.  The  computational  grid  must  be  fine 
enough  to  resolve  the  rotational 
structures. 

Feature  6  implies  a  direct  numerical 
simulation  of  the  conservation  equations. 


without  modeling  the  fluctuating  flow 
components.  Fortunately,  feature  7  can  be 
satisfied  because  these  structures  are  large 
and  can  be  resolved  on  computational  grids 
that  are  orders  of  magnitude  coarser  than 
those  needed  to  resolve  molecular  transport 
effects. 

Shear  layers  occurring  in  most  practical 
applications  will  experience  three-dimen¬ 
sional  perturbations,  and  eventually  the 
fluctuating  kinetic  energy  will  be  shared 
over  all  of  these  dimensions.  Hence,  the 
principal  limitation  of  the  present  calcula¬ 
tions  is  the  two-dimensional  flow  approxi¬ 
mation. 

Nevertheless,  the  present  large-Reynolds- 
number  approach  (i.e.,  the  direct  solution 
of  the  conservation  laws  of  gas  dynamics) 
offers  one  technique  for  extending  the 
numerical  simulations  of  fluctuating  shear 
layers  to  the  nonsteady,  compressible  flow 
regime.  Such  numerical  solutions  are 
valuable  for  at  least  three  reasons:  ( 1 )  they 
provide  idealized  solutions  against  which  to 
compare  real  flows:  (2)  the  problem 
parameters  (e.g..  initial  and  boundary 
conditions)  can  be  controlled  more  easily 
than  in  experiments:  (3)  the  flow  field  can 
be  sampled  nonintrusively  and  easily 
analyzed.  Such  comparisons  between  ideal¬ 
ized  calculations  and  experimental  data 
cannot  help  but  increase  our  understanding 
of  the  fluid  dynamics  of  such  complex 
flows. 

In  closing,  it  should  be  pointed  out  that  the 
organized  structures  of  the  shear  layer 
flows  presented  here  can  be  viewed  as 
noteworthy  examples  of  self-organized 
response  of  nonequilibrium  systems  as 
described  by  Nicolis  and  Prigogine  (1977) 
and  Prigogine  (1980).  Whereas  low- 
Reynolds-number  shear  layers  smooth  out 
tangential  velocity  discontinuities  by  molec¬ 
ular  diffusion,  large-Reynolds-number 
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shear  layers  remove  the  velocity  nonequili¬ 
brium  by  forming  large-scale  rotational 
structures.  The  initial  and  boundary 
conditions  of  the  problem  (i.e.,  shear  flow) 
place  the  system  in  a  state  corresponding 


to  the  nonlinear  range  of  nonequilibrium 
dynamics.  The  system  is  unstable,  and  the 
flow  evolves  to  a  new.  more  stable  state 
which  is  considerably  more  complex  yet 
self-organized. 
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